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I. ASPECTS OF LITHOSPHERIC EVOLUTION ON VENUS 
II. THERMAL AND COLLISIONAL HISTORIES OF CHONDRITE PARENT BODIES 

by 

Robert E. Grimm 

Submitted to the Department of Earth, Atmospheric, and Planetary Sciences on 
October 25, 1988, in partial fulfillment of the requirements of the Degree of 
Doctor of Philosophy in Geophysics 


ABSTRACT 

This thesis consists of two principal sections which address the geological evo- 
lution of distinctly different kinds of solar system objects. Venus, the second largest 
of the terrestrial planets, has been observed over the past decade by orbital radars 
on both American and Soviet spacecraft. These surface measurements provide 
clues to the structure and evolution of the lithosphere. The parent bodies of chon- 
dritic meteorites, thought to resemble asteroids, represent the other end of the size 
spectrum of terrestrial objects. Their early thermal and collisional histories may be 
constrained by the chemical and textural record preserved in meteorite samples. 

Impact craters on Venus have been observed by the Soviet Venera 15/16 
spacecraft. We present a formalism by which the size-frequency distribution of 
impact craters may be used to estimate upper bounds on the mean global rates of 
volcanic resurfacing and lithospheric recycling on that planet over the past several 
hundred million years. The principal assumptions are that an upper bound on the 
rate of crater production on Venus may be estimated, that craters are volcanically 
obliterated only when flows completely cover their rims, and that the mean rate of 
lithospheric recycling is inversely proportional to the surface age. The impact 
crater density reported from Venera 15/16 observations, if valid for the entire Venus 
surface, then indicates a mean volcanic flux no greater than 2 km 3 /y, correspon- 
ding to a maximum average rate of resurfacing of about 4 km/b.y. The fraction of 
global heat loss due to such a rate of volcanic resurfacing is negligible. For the 
lowest estimated mean crater retention age of the surface of Venus imaged by 
Venera 15/16, the rate of lithospheric recycling on Venus does not exceed 1.5 
km 2 /y, corresponding to 25% of the global heat loss. More likely estimates of the 
mean surface age limit the contribution to global heat loss of any lithospheric 
recycling to be less than 10%. These results support the hypothesis that simple 
conduction dominates heat transport at lithospheric levels in the Venus interior. 

Changes in the morphology of impact craters with time may also be used to 
investigate the interior. Because of the high surface temperature on Venus and the 
strong temperature dependence of strain rate in silicates, solid state creep may be 
an important mechanism for the reduction of topographic relief on that planet. On 
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the basis of models for the viscous relaxation of impact crater topography we derive 
constraints on the thickness of the crust and the mean lithospheric thermal gradient 
beneath craters observed by Venera 15/16. From the mean and variance of ob- 
seved crater depth as a function of crater diameter on Venus and from estimates of 
the initial depths of fresh impact craters on Venus obtained by scaling from lunar 
observations, we formulate a statistical test for the maximum amount of viscous 
relaxation that is consistent with the observations at a specified confidence level. 
We develop a general formulation for gravity-driven flow in a linearly viscous fluid, 
incorporating the densities and temperature-dependent effective viscosities of dis- 
tinct crust and mantle layers. The statistical comparison of predicted and observed 
depths yields linked upper limits to crustal thickness H and thermal gradient dT/dz. 
Since dT/dz can be estimated from global heat loss arguments or thermal models, 
an upper bound to H may be derived. The preserved relief of the largest craters 
constrains H to be less than or equal to 10-20 km. Because the craters with meas- 
ured depths appear to be representative of a larger population of impact structures 
in terms of regional elevation, geological unit, and degree of preservation of mor- 
phological detail, extrapolation of this result to lowlands and rolling plains regions 
on a global basis is not unreasonable. Such an extrapolation, together with iso- 
static considerations, yields an upper limit to the crustal volume of Venus of 10 10 
km 3 . This value is an order of magnitude less than the time-integrated volume of 
crust produced on Earth and implies either that the average rate of crustal gener- 
ation has been much smaller on Venus than on Earth or that some form of crustal 
recycling has occurred on Venus. 

Aphrodite Terra, the largest highland region of Venus, is a likely site of mantle 
upwelling and active volcanism and extensional tectonics. Crumpler and Head 
have proposed that Aphrodite is a divergent plate boundary and cite as support of 
this model an organized system of lineaments held to be analogs to oceanic 
fracture zones; small-scale bilateral symmetry of topographic elements, hypothe- 
sized to be rifted and separated relief analogous to some terrestrial oceanic 
plateaus; and subsidence of topography proportional to the square root of distance, 
consistent with that expected for a divergent thermal boundary layer. We under- 
take quantitative tests of each of these assertions, and we compare the results to 
similar tests of the Mid-Atlantic Ridge. We find that, apart from the long-wavelength 
symmetry of broadly elevated regions, there is no evidence for regional bilateral 
symmetry of features several hundred kilometers in size on either planet. The fit of 
the topography of Aphrodite to a that of a thermal boundary layer is in general 
much poorer than for the Earth, and so other mechanisms such as dynamic uplift or 
crustal thickness variations must dominate the topography. The broad saddle- 
shaped region between Thetis and Atla Regiones shows the best fit to the root- 
distance relation, yielding apparent spreading half-rates of a few centimeters per 
year, but the statistics are quite sensitive to the distance range analyzed. Calcula- 
tion of a single pole of relative motion for the entire postulated system of transform 
faults shows that the inferred fracture zone traces are not consistent with a simple 
two-plate model, regardless of past pole motions. A multiple-plate geometry is 
required, and one or more of the lineaments must be a plate boundary if the plate- 
divergence hypothesis is correct. Such a boundary would be discernable by 
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geological evidence for non-transform motion in forthcoming Magellan radar 
images. If, however, in the absence of such evidence, the lineaments themselves 
are nevertheless verified, then rigid-plate rotations are invalid, and a nonrigid 
model must be adopted in which deformation is broadly accommodated and the 
lineaments are surface traces of mantle convective flow. A model dominated by 
vertical tectonics, wherein horizontal mantle divergence is largely decoupled from 
the surface, is possible only if the lineaments are unconfirmed by Magellan . 

Ordinary chondrite meteorites show textural and chemical patterns indicative of 
varying intensities of thermal metamorphism. The conventional "onion-shell" 
model, which envisions highly metamorphosed material in the core and less 
intensely heated rocks near the surface, predicts an inverse relation between peak 
temperature and cooling rate, but none has been observed. Scott and Rajan 
devised a "metamorphosed-planetesimal" model to explain this discrepancy, 
whereby heating occurs in planetesimals a few kilometers in radius which then 
accrete to form 100-km-radius parent bodies. Cooling rates are then randomly 
controlled by burial depth. Thermal and collisions constraints on the metamor- 
phosed-planetesimal are examined here, and the model is found to be applicable 
only to highly insulating, 26 AI-rich planetesimals that remain closely aggregated 
upon accretion. An alternative model is presented here, in which onion-shell 
parent bodies are collisionally fragmented during metamorphism and then gravi- 
tationally reassembled. If reassembly times are short, then cooling rates would be 
determined by burial depth in the reaccreted parent body. This model, unlike 
previous ones, can explain both coherent and incoherent cooling of breccia clasts, 
by collisions during or after metamorphism, respectively. 

Carbonaceous chondrites have, in general, not experienced the thermal 
metamorphism seen in ordinary chondrites but have been aqueously altered within 
their parent bodies. From chemical and textural data on these meteorites, and from 
studies of collision mechanics, we pose two hypotheses for the aqueous alteration 
environment. In the first model, alteration occurs uniformly throughout the parent 
body interior; in the second, alteration occurs in a post-accretional surface regolith. 
Both models are based on the assumptions of an initially homogeneous mixture of 
ice and rock and heating by decay of 26 AI. Under the interior-alteration model, 
linked bounds on the initial ice-to-rock ratio and 26 AJ abundance may be found that 
satisfy peak temperatures derived from oxygen isotope studies. Additional 
constraints imposed by theJnferred water volume consumed by the alteration 
reaction and the total water volume that exchanged oxygen isotopes with host 
rocks are best explained if alteration occurred in a regolith. We show quantitatively 
how liquid water may be introduced there by hydrothermal circulation, by diffusion 
of vapor from below, or by venting due to fracture when interior pore pressures 
exceed the parent body strength. Retention of primordial ice is probably not limited 
by sublimation or by collisional comminution, but by shock vaporization. If large, 
C-type asteroids are representative of carbonaceous chondrite parent bodies, they 
may contain significant quantities of ice. 

Thesis Supervisor: Sean C. Solomon, Professor of Geophysics 
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CHAPTER 1 : 


INTRODUCTION 



The solar system is a motley family. Its members vary vastly in size, from Jupiter, 
the glaring master of the night sky, to the countless dust particles which comprise 
the hazy zodiacal light. Between these extremes are those solid bodies which have 
undergone distinct structural and chemical changes with time, that is, they have 
geological histories. These objects range in size from Earth and Venus, with dia- 
meters of several thousand kilometers, to meteorite parent bodies, whose diameters 
are inferred to have been several hundred kilometers or less. This thesis inves- 
tigates several aspects of the physical evolution of objects at both boundaries of 
geological inquiry, the planet Venus and the parent bodies of chondritic meteorites. 

The surface geological histories of terrestrial objects are determined principally 
by the competition between endogenic processes, manifested by volcanism and 
tectonism, and the exogenic imprint of impact cratering (erosion plays a secondary 
role). The earliest geological evolution of terrestrial objects was dominated by 
collisions during their accretion. Continued impact cratering, with decreasing 
intensity over solar system history, has affected all planetary surfaces. The ubiquity 
of impact craters makes them not only the products of an easily studied geological 
process in themselves, but variations in their areal density are excellent measures 
of the relative ages of different regions of a planetary surface, although absolute 
ages are more difficult to assign. Moreover, detailed examination of the size- 
frequency distribution and morphology of impact craters can be used to infer 
mechanisms and rates of endogenic modification, which can in turn yield clues on 
the state of the interior. 

Heating of a planetary interior drives volcanism and tectonism at the surface. 
For the planets and larger moons, the main heat source is radioactive decay. 
Heating is proportional to the planetary volume for such an internal heat source. 
However, radiation is proportional to the surface area, which leads to the so-called 
"square-cube law": since the volume-to-area ratio increases linearly with planetary 
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size, larger planets will retain heat better and have more vigorous and longer- 
lasting geologic histories. This simple and perhaps well-worn rule nevertheless 
explains the approximate correlation between size and surface age inferred from 
impact crater density: the Moon and Mercury, heavily cratered, have been largely 
inactive for billions of years; Mars, partially cratered but marked with enormous 
shield volcanoes and deep rifts, attests to activity within the past few billion years; 
Earth, with only a smattering of craters largely on old continental shields, remains 
strongly volcanically and tectonically active today. Using size as a guide to thermal 
longevity leads to predictions that Venus should be presently volcanically and tec- 
tonically active, whereas small meteorite parent bodies should be long dead. The 
active early histories of these objects are nonetheless preserved in the chemistry 
and texture of meteorites. Understanding the earliest geological processes of the 
solar system and the different geological histories of the two largest terrestrial 
objects are important goals in planetary science. 

The relationship between impact and thermal histories is explored along several 
avenues in this thesis. In the first two investigations, impact craters are used as 
tools to probe the outer tens of kilometers of Venus. Chapter Two combines the 
size-frequency distribution of impact craters on Venus with estimates of the recent 
cratering rate to calculate the maximum rates at which impact craters can be obliter- 
ated by volcanism or lithospheric recycling. The heat loss by each mechanism may 
be derived from the crater obliteration rate, therefore leading to identification of the 
dominant mode of heat loss at lithospheric levels on Venus. 

Chapter Three analyzes the endogenic modification of impact craters on Venus 
by viscous relaxation and derives constraints on the rheological and thermal 
structure of the interior imposed by observed crater morphology. Viscous relax- 
ation, a process analogous to terrestrial postglacial rebound wherein gravity-driven 
flow gradually eradicates topographic relief, is expected to more important on 
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Venus because of the high surface temperature. The viscous relaxation model 
developed in this thesis depends principally on the local crustal thickness and 
thermal gradient. Statistical comparison of the observed depths of craters against 
model predictions of the depth of partially relaxed craters is used to establish linked 
bounds on these parameters. 

The tectonic styles resulting from specified modes of heat transport are investi- 
gated in Chapter 4 for Aphrodite Terra, the largest highland region on Venus. Three 
alternative models for the interaction of mantle convection with the surface are 
considered (vertical tectonics, rigid-plate divergence, and distributed deformation). 
These models differ in the strength of coupling of interior flow to the surface and the 
response of the surface to these stresses. A rigid-plate model is the most straight- 
forward to test, as the associated topography and tectonics have been well studied 
for terrestrial ocean basins. Departures from the predictions of this model may be 
used to assess the alternative hypotheses. Further predictions necessary to 
discriminate between these models may be resolved by forthcoming high-resolution 
radar imaging by the Magellan spacecraft. 

Impact and thermal processes are again linked in Chapter 5, which examines 
the metamorphic and collisional history of ordinary chondrite parent bodies. 
Chemical and textural properties of meteorites yield information both on the peak 
temperatures during metamorphism and subsequent cooling rates. Because no 
simple relation has been observed between these parameters for ordinary chon- 
drites, some workers have suggested that metamorphic signatures were fixed in 
small planetesimals, and that cooling was randomly controlled by the timing and 
location of accretion into the final parent bodies. Some implicit thermal and col- 
lisional assumptions in this model are clarified here. These considerations also 
suggest a new model in which parent bodies are collisionally fragmented and 
gravitationally reassembled during the metamorphic interval. 



The thermal history of carbonaceous chondrite parent bodies has heretofore 
been largely ignored because these primitive meteorites have not been thermally 
metamorphosed. However, aqueous alteration has profoundly affected their 
chemistry and texture, and some heat source must be responsible for these effects. 
In Chapter 6, thermal models for carbonaceous chondrites are tested under the 
alternative assumptions that aqueous alteration occurred in the interiors of parent 
bodies or within surficial regoliths. The principal constraints on these models derive 
from oxygen isotope systematics. The thermodynamic and transport properties of 
H 2 0 bear strongly on the thermal evolution of interior-alteration models. Several 
specific mechanisms for introducing water to the regolith are suggested and quanti- 
tatively supported. The secular retention of ice on asteroids is also discussed, and 
the implications for exploration are assessed. 

Chapters Two and Three were cowritten with S.C. Solomon and appeared in 
Geophys. Res. Lett., 14, 538-541, 1987 and J. Geophys. Res., 93, 11911-11929, 
1988, respectively. Chapter Four, also written with Solomon, will be submitted to 
the Journal of Geophysical Research. Chapter Five was published in J. Geophys. 
Res., 90, 2022- 2028, 1 985. Chapter Six, cowritten with H.Y. McSween, will be 
submitted to Icarus. 
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CHAPTER 2: 

LIMITS ON MODES OF LITHOSPHERIC HEAT TRANSPORT 
ON VENUS FROM IMPACT CRATER DENSITY 
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INTRODUCTION 

Characterization of the dominant mode of lithospheric heat transport on Venus 
remains an important step in understanding the internal evolution of the planet 
[Solomon and Head, 1 982]. Among the other solid planets and satellites of the 
solar system are examples of bodies on which lithospheric heat transport is 
dominated by conduction (Moon, Mercury, Mars), volcanism (lo), and lithospheric 
recycling (Earth). Figure 2.1 is a schematic ternary diagram showing the estimated 
relative contribution of each of these heat loss mechanisms for these bodies. 
Images of the Venus surface from Earth-based radar [Campbell et ai, 1983, 1984] 
and Venera 15/16 [Barsukov et ai, 1986; Basllevsky et al., 1986] have revealed 
evidence for widespread volcanism and tectonism but also a significant number of 
impact craters [Campbell and Bums, 1980; Ivanov et ai, 1986]. On the basis of 
simple models of impact crater production and obliteration, we present a formalism 
to estimate upper bounds on the contributions to lithospheric heat transport by 
volcanism and lithospheric recycling from the observed density of impact craters on 
the Venus surface. 

CRATERING 

Impact crater densities and crater retention ages are usually expressed relative 
to a reference surface, typically the lunar maria [BVSP, 1981]. Under the 
assumption that the same population of planetesimals has been responsible for 
cratering throughout the inner solar system, surfaces of the same age on the Moon 
and Venus should have crater density distributions differing only by a multiplicative 
scaling factor derivable in principle from the relative impact probabilities and 
velocities of planet-crossing asteroids and comets [ Hartmann , 1977; BVSP, 1981]. 
Of course, the atmosphere of Venus will act to reduce the production of impact 
craters, particularly at smaller diameters, because of the deceleration and break-up 
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of incident meteors [Tauber and Kirk, 1976; Kahn, 1982]. Assuming that the 
cratering rate on each body has been constant since emplacement of the lunar 
reference surface, the expected number of craters per area produced in the 
diameter increment 0 to V2D per time on an atmosphere-free Venus can be 
expressed by 


F V (D) = (f N rT /t m )D-b (2.1) 

where f is the ratio of the crater production rate on Venus to that on the Moon, t m 
and N m are the age and incremental crater density at unit diameter, respectively, of 
the lunar reference surface, and b is a constant. Similarly, the average crater 
retention age t v of the Venus surface is [Hartmann, 1977] 


t v * t m Ny/(f Nm) 


( 2 . 2 ) 


where N v is the incremental crater density at unit diameter for Venus. 

Given the uncertainties in the lunar cratering record and particularly in the 
scaling of crater production from the Moon to Venus, it is appropriate to consider a 
range of plausible cratering rates on Venus. Several rates are presented in Table 
2.1. For each estimate, the lunar mare regions surrounding the Apollo 12 and 15 
landing sites are taken as the reference surface. Mare basalts from Apollo 12 and 
15 are among the youngest dated lunar rocks, and the assumption of a near- 
constant cratering flux since the emplacement of these mare units is a reasonable 
one [BVSP, 1981]. Uncertainties in both crater density and mare emplacement age 
[BVSP, 1981] are represented. We consider two independent estimates, with 
associated uncertainties, in the crater production scaling factor f. The first estimate, 
f m 0.8 - 2.0 [Hartmann, 1977; BVSP, 1981], has been applied to the cratering 
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record of Venus by a number of other workers [Masursky et at., 1980; Campbell 
and Burns, 1980; Ivanov et ai, 1986]. The second and more recent estimate, f = 
1.2 - 4.2, is derived from Schaber et al. [1987], who recalculated impact 
probabilities and concluded that impact rates on Earth and Venus are nearly 
identical. In particular, they considered that the estimate of Grieve [1984] for the 
cumulative production function for craters with D > 20 km during the past 120 m.y. 
on Earth, (5.4 ± 2.7)x10 -15 knrV 1 , is valid for Venus. We compare this value to the 
cumulative cratering rate for D > 4 km at the Apollo 12 and 15 sites, correcting for 
the difference in cutoff diameter using b = 1.8 [BVSP, 1981]. The equivalent f is 
then the ratio of crater production on the Earth to that on the Moon. For both sets of 
estimates of f, we take the unit diameter to be 1 km. This choice is largely for ease 
of evaluation of equation (2.1) when D is expressed in kilometers; we recognize 
that N m and N v have meaning only as values extrapolated from crater densities at 
diameters sufficiently large so that the effects of secondary craters on the Moon and 
of atmospheric shielding on Venus are unimportant. The range in incremental 
crater production rate on Venus so extrapolated to D = 1 km, F v (1) = f is then 
as given in Table 2.1. 

A mean crater retention age for Venus can be calculated from equation (2.2) for 
a given density of impact craters on the Venus surface. Different assumptions have 
been made in the literature as to which surface features visible in radar images are 
likely to be of impact origin. Impact craters have been identified by Campbell and 
Burns [1980] in Earth-based radar images as circular areas of low backscatter 
surrounded by an annular region of high contrast terrain; such features sometimes 
display a central bright spot. Masursky et ai [1980] drew attention to a number of 
large, circular, radar-dark depressions and suggested they may also be impact 
features. The 1-4 km resolution of Venera 15/16 images allows more confident 
identification of such classic characteristics of impact craters as elevated rims, 
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central peaks, and ejecta blankets [Ivanov et al, 1986], although there remain 
numerous approximately circular features of uncertain origin [Barsukov et al.. 
1986]. The Venera 15/16 images thus provide the best current estimate of crater 
density statistics for the northern quarter of the planet [Ivanov et al., 1986]. We 
note, however, that the crater density distribution reported by Campbell and Bums 
[1980] for a different but overlapping region is similar for craters in the diameter 
range 80-140 km. 

We have estimated N v from a linear regression on the Venera 15/16 
incremental crater density data (log-log plot) of Ivanov et al. [1986] for D > 23 km, 
with b fixed at 1 .8. This estimate, with associated uncertainties, can be substituted 
into equation (2.2) to define a range of possible crater retention ages (Table 2.1). A 
crater retention age of about 1 b.y. is obtained if f si , in agreement with the mean 
ages reported by Campbell and Burns [1980] and Ivanov etal. [1986] under similar 
assumptions. When all uncertainties are considered, however, the range in 
possible mean surface ages is quite large, from about 100 m.y. to over 2 b.y. 

VOLCANIC RESURFACING 

Impact craters, of course, can be obliterated by many processes, including 
erosion, sedimentation, viscous relaxation, and tectonic disruption, as well as 
volcanic resur- facing and lithospheric recycling. Several studies [Garvin et al., 
1984; Head et al., 1985; Barsukov et al., 1986] suggest that the effects of erosion 
and sedimentation on Venus are unimportant at the scale of features considered 
here. As a working hypothesis, we first assume that volcanism is the major 
resurfacing process and is equally likely anywhere on the Venus surface. This 
hypothesis is at least partly supported by the widespread occurrence of smooth 
plains, flow features, and volcanic cones, domes, and calderas [Campbell et al., 
1984; Barsukov et al., 1986]. The combined effects of other obliteration processes 



20 


serve only to decrease the contribution of volcanism required, which ensures that 
our estimate of volcanic flux will be an upper bound. 

To maintain a conservative approach, we assume that a crater must be buried 
by volcanic flows to the top of its rim in order that it not be observable in a radar 
image. We estimate the rim height h of fresh craters on Venus by applying a gravity 
scaling factor to the power-law relation between h and diameter D for lunar craters 
[Pike, 1977]. Two different scaling rules are possible, following suggestions that 
either the depth (and, by inference, the rim height) of complex craters or the 
transition diameter between simple and complex crater morphology varies 
inversely with gravity [Grieve and Robertson, 1979; Pike, 1980]. We choose the 
latter scaling rule here because it predicts higher rims and hence requires a 
greater amount of volcanic resurfacing to achieve crater obliteration. The 
prediction for Venus rim heights, in km, is then h = 0.085 D 0 - 4 . 

For a specified time t since initiation of resurfacing at a constant rate of volcanic 
burial ds/dt, there exists a critical crater diameter D c at which the rim height is equal 
to the total thickness of volcanic flows: h(D c ) = t (ds/dt). Craters with D > D c are 
unaf-fected by the resurfacing process and reflect the crater density prior to 
resurfacing; craters with D < D c are continually obliterated and replaced. The 
crater density in the latter state is in equilibrium and depends on the ratio of the rate 
of cratering to that of obliteration [Jones, 1974] 

N(D) = F(D) h(D)/(ds/dt) (2.3) 

where N is the number of preserved craters per area in the diameter interval D to 
V2D. The volcanic resurfacing model is schematically illustrated in Figure 2.2. For 
a specified crater density, a larger volcanic flux is allowed for a larger crater 
production rate. For F(1) = l^xIO-^knrrfy* 1 . the largest value in Table 2.1, the 



21 


predicted equilibrium crater densities for steady volcanic fluxes of 1 , 2, and 4 km 3 /y 
are as shown in Figure 2.3. If we take the diameter of the largest reported crater in 
the Venera images (140 km) as the critical diameter D c> the minimum periods of 
resurfacing so that equation (2.3) is valid over the diameter range spanned by the 
Venera data are 280, 140, and 70 m.y., respectively. 

On the basis of Figure 2.3, an upper limit to the mean volcanic flux consistent 
with the Venera observations, if representative of the entire planet, is 2 km 3 /y. This 
is equivalent to an average volcanic resurfacing rate ds/dt of 4 km/b.y. Of course, 
volcanic extrusion has not been steady but has likely fluctuated about a gradually 
decreasing function of time. The inferred upper bound on the rate of resurfacing on 
Venus remains valid, however, as such a volcanic history would have eliminated 
more craters than in a constant-flux model. 

An independent geochemical test of this result has been obtained by Fegley 
and Prinn [1988], who concluded that a volcanic eruption rate of 0.4 to 1 1 km 3 /y is 
required to replenish atmospheric S0 2 lost to chemical weathering of surface 
rocks. Their best estimate, using surface compositions reported by Venera and 
Vega landers, is 1 km 3 /y, quite consistent with the value reported here. 

The rate of heat Q v delivered to the surface by a volcanic flux dV/dt is 

Q v = p(CpAT + H f ) (dV/dt) (2.4) 

where C p is the heat capacity, p is the density, Hf is the heat of fusion, and AT is the 
temperature drop across the lithosphere. Using AT = 1015 K [Kaula and Phillips, 
1 981 ], p » 2.8 g/cm 3 , C p = 1 .2 J/g K, and Hf = 400 J/g, the rate of heat loss due to a 
volcanic flux of 2 km 3 /y is 2.9 x 10 11 W, or less than 1% of the nominal rate of 
global heat loss of 3.4 x 10 13 W [Solomon and Head, 1982]. This rate of heat 
transport is negligible on a planetary scale, though such a result does not preclude 
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the possibility of a significant component of heat transport by shallow conduction 
above igneous intrusions [Head and Wilson, 1986]. Volcanism on Earth, 
principally associated with lithospheric divergence, delivers an order of magnitude 
more heat per time than inferred here for Venus. However, the terrestrial volcanic 
contribution to heat loss is still nearly an order of magnitude smaller than that lost 
due to creation and cooling of the lithosphere [Sclater et at., 1980], which 
distinguishes plate recycling as a separate mode of heat loss (Figure 2.1). 

LITHOSPHERIC RECYCLING 

We next consider as an alternative working hypothesis that the principal 
mechanism of impact crater obliteration is tectonic disruption due to lithospheric 
recycling. Such recycling may be considered to occur either by foundering or by 
subduction, though we recognize that evidence for neither process has yet been 
observed in images of the Venus surface [Barsukov et ai, 1986; Basilevsky et ai, 
1986]. If the regions imaged by Venera 15/16 are globally representative, then the 
lowest estimate of the mean crater retention age (130 m.y.) in Table 2.1 can be 
converted into an upper bound on the rate of lithospheric recycling. If we assume 
that there are no preferentially preserved "islands" of more ancient terrain, the rate 
of lithospheric recycling is inversely proportional to the mean age of the surface. 
The mean age of oceanic lithosphere on Earth is approximately 60 m.y. [Sclater et 
ai, 1980], which would imply under the above assumptions that any lithospheric 
recycling on Venus is at least a factor of 2 slower than on Earth. The terrestrial rate 
of creation of new lithosphere is 3 km 2 /y [Parsons, 1981]; the maximum rate on 
Venus would thus be 1.5 km 2 /y. The presence of ancient terrain uninvolved in the 
recycling process would seriously bias this estimate; for instance, the mean age of 
the Earth's surface including the continents is about 400 m.y. [Sclater et ai, 1980]. 
However, no systematic large scale variations in crater density are evident in the 
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Venera images [Barsukov et at., 1986]. In particular, no regions of heavily cratered 
terrain have been identified in the portions of the Venus surface well imaged by 
Arecibo and Venera. 

The rate of heat loss due to lithospheric recycling is approximately 

Q = 0.5 p C p L AT (dA/dt) (2.5) 

where L is the thickness of the thermal lithosphere and dA/dt is the areal rate of 
creation of new lithosphere. With dA/dt = 1 .5 km 2 /y, C p and AT as above, p = 3.3 
g/cm 3 and L = 94 km [Kaula and Phillips, 1981], Q = 9 x 10 12 W, or 25% of the 
nominal Venus heat loss. The fractional heat loss corresponding to the more likely 
estimates of cratering rate and mean crater retention age in Table 2.1 are only 
3-9%, a small fraction of the Venus heat budget. These estimates of the rate of 
Venus heat loss contributed by lithospheric recycling are comparable to the upper 
bound (15%) derived by Kaula and Phillips [1981] using different arguments. It is 
evident that the situation on Venus is in strong contrast to the Earth, where plate 
recycling accounts for approximately 65% of the heat loss from the interior [Sclater 
etal., 1980]. 

CONCLUSIONS 

The density of impact craters in Venera 15/16 radar images of the Venus 
surface, if representative of the entire planet, limits the average rate of volcanic 
resurfacing to be less than 2 km 3 /y, corresponding to less than 1% of the global 
heat loss. On similar grounds the rate of lithospheric recycling on Venus is less 
than 1 .5 km 2 /y and probably less than 0.5 km 2 /y; these figures correspond to 25% 
and 9% of the global heat loss, respectively. The position of Venus on the ternary 
diagram of heat loss mechanisms as estimated by this work is given in Figure 2.4. 
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These results suggest that heat loss at lithospheric levels in Venus is dominated by 
conduction. Substantial lateral variation in vertical thermal gradients, associated 
with lithospheric thinning [Morgan and Phillips, 1983] or magmatic intrusions [Head 
and Wilson, 1986], may nonetheless be present. Detailed characterization of 
volcanic, tectonic, and impact processes on Venus must await the improved 
resolution and coverage of radar images to be obtained during the Magellan 
mission and will provide tests of the simple models and arguments presented here. 

APPENDIX 2.A: Effect of geographic variations in crater density on the 
estimated rate of lithospheric recycling. 

Both of the approaches discussed above require that there be no significant 
heterogeneity in crater density with geographic location. On the basis of a x 2 test 
of distribution, Basilevsky et al. [1987] found no reason to reject the null hypothesis 
that impact craters are randomly distributed over the area imaged by Venera 15/16. 
However, Head et al. [1988] concluded that crater density over this portion of the 
Venus surface increases as a function of latitude. They took this variation as 
supportive of a model of equatorial crustal spreading [Head and Crumpler, 1988; 
see Chapter 4 of this thesis]. We have confirmed this latitudinal trend in the Venera 
data [Basilevsky et al., 1987] with a x 2 test, although craters identified from the 
Arecibo data [Campbell and Burns, 1980] in a different region of Venus do not 
support a similar result. 

If such a zonal variation in surface age is present and is due to creation of 
young surface near the equator, its effect on the mean surface age and contribution 
to heat loss by lithospheric recycling may be estimated. Under the simplest 
assumption that impact crater density increases linearly with latitude from a value 
of zero near the equator, the crater density N over a latitudinal band 0! to 0 2 is 
proportional to 
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J 0 cos 0 d0 

N 

e 2 

J cos 0 d0 

e i 

From equation (2.6), we find that the mean crater density over the portion of the 
surface imaged by the Venera orbiters (approximately 30°-90°N ) would be 1 .6 
times larger than that of the entire northern hemisphere if crater density increases 
linearly with latitude. Assuming a symmetrical crater distribution for the southern 
hemisphere leads to a corresponding upper limit to the fractional heat loss of 
-35%. More likely estimates of the cratering rate given in Table 2.1 still limit this 
contribution to less than 1 5%, however. 
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TABLE 2.1. Incremental crater production function F V (D) at D = 1 km 
and implied crater retention age t v for Venus 


f 

109 y 

N m> 

1 0* 4 km* 2 

F v (1). 

10-13 km* 2 y 1 

N v , 

1 0 -4 km* 2 

t 

109 y 

Scaling factor f from Hartmann [1 977] and BVSP [1 981 ] 


0.8 a 

3.4 b 

4.4a 

l.ia 

2.4 b 

2.3 b 

1.0 

3.2 

7.2 

2.2 

2.1 

0.94 

2.0 b 

3.1 a 

11 P 

6.8 b 

1.9a 

0.27a 

Scaling factor f from Grieve [1 984] and Schaber et al. [1 987] 


1.4 a 

3.4 b 

4 . 4 a 

1.9a 

2.4 b 

1.3 b 

2.8 

3.2 

7.2 

6.3 

2.1 

0.34 

4.2 b 

3.1® 

11 b 

14 b 

1.9a 

0.13a 


a Approximate lower bound for this quantity. 
b Approximate upper bound for this quantity. 
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Figure 2.1. Schematic ternary diagram showing the relative contributions of each 
mechanism of lithospheric heat transfer toward global heat loss on 
solid planets and satellites, after Solomon and Head [1982]. For the 
Moon, Mercury, and Mars, thermal conduction carries all of the 
lithospheric heat flow. The position of Earth, on which the majority of 
lithospheric heat transport occurs by plate recycling, does not include 
the contribution to global heat loss from continental crustal radioactivity. 
Hot-spot volcanism dominated heat loss on lo, but the relative fraction 
contributed by conduction is not well known; at the position shown, 
the conductive heat flow on lo is equal to the Earth's average heat flux. 
The placement of Venus on this diagram is uncertain. 

Figure 2.2. Schematic view of how impact crater density limits the rate of volcanic 
resurfacing. The crater density N(D) is in equilibrium for D > D c (see 
text). The crater-production function F(D) is the expected number of 
craters within a given diameter increment produced per area per time. 
For assumed values of F(D) and h(D), the observed crater distribution 
limits ds/dt. 

Figure 2.3. Incremental crater density versus diameter obtained from Venera 1 5/1 6 
observations [Ivanov et ai, 1986] compared with equilibrium crater 
densities for constant assumed global volcanic fluxes of 1 , 2, and 4 
km 3 /y (solid lines). 

Figure 2.4. Position of Venus on the ternary diagram for mechanisms of 
lithospheric heat transport, as constrained by impact crater density on 
the portion of the Venus surface imaged by Venera 15/16. 
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CHAPTER 3: 


VISCOUS RELAXATION OF IMPACT CRATER 
RELIEF ON VENUS: CONSTRAINTS ON CRUSTAL 
THICKNESS AND THERMAL GRADIENT 
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INTRODUCTION 

Radar images of Venus have revealed a planet with a variety of tectonic and 
volcanic features [Campbell et al., 1 983, 1 984; Barsukov et a/., 1 986]. The Venus 
surface has also preserved a number of large impact craters [Campbell and Burns, 
1980; Ivanov et al., 1986; Basilevsky et al., 1987]. The density and distribution of 
these craters provide crucial information for understanding the geological evolution 
of the planet [Ivanov et al., 1986; Schaber et al., 1987a; Grimm and Solomon, 
1987]. The endogenic modification of impact craters can also yield key constraints 
on internal planetary processes. One such modification mechanism is the viscous 
relaxation of crater topography. This process is thought to have been important for 
the Moon [Danes, 1965; Hall et a!., 1981; Solomon et al., 1982a; Baldwin, 1987] 
and the Galilean satellites [Parmentier and Head, 1981 ; Passey and Shoemaker, 
1982 \ Thomas and Schubert, 1986, 1987]. On Venus, considerable viscous relax- 
ation of topography has been anticipated because of the high surface temperature 
and the strong temperature dependence of strain rate [Weertman, 1979; Solomon 
etal., 1982b, Stephens etal., 1983]. Venera 15/16 radar altimetric profiles across 
impact craters 30 - 140 km in diameter, however, show crater depths in excess of 
several hundred meters, suggesting that the extent of relaxation of relief for these 
features has been small [Ivanov et al., 1986; Basilevsky et al., 1987]. 

In this paper we employ the observed depths of impact craters to test models for 
the viscous relaxation of surface topography on Venus. We utilize a general 
numerical formulation of gravity-driven flow in a linearly viscous fluid, and we 
incorporate distinct densities and temperature-dependent flow laws for crustal and 
mantle material. Given an estimate for the initial topography of impact craters on 
Venus, we compare the observed topography with that predicted by viscous 
relaxation for a range of possible values for crater age, crater diameter, crustal 
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thickness, and mean lithospheric thermal gradient. On the basis of a statistical 
approach, we determine confidence levels for rejecting combinations of these 
parameters. These calculations suggest that the mean thickness of the crust over a 
substantial fraction of the planet is probably about 10-20 km. We consider several 
possible implications of this result for the evolution of Venus. 

IMPACT CRATERS ON VENUS 

Impact craters on Venus were first identified in Earth-based radar images of the 
surface as circular areas of low backscatter surrounded by an annular region of 
high reflectivity; the circular radar-dark regions sometimes contained a central 
bright spot [Campbell and Burns, 1980]. The orbital radar images obtained by 
Venera 15 and 16, with a resolution of 1-2 km, currently allow the best identification 
of impact craters over the imaged portion of the planet [Ivanov et ai, 1986; 
Basilevsky et ai, 1987]. Distinct characteristics diagnostic of impact craters, such 
as elevated rims, central peaks, and ejecta blankets, can be resolved for many 
craters, although there are also many approximately circular features of uncertain 
origin [Barsukov et ai, 1986]. A total of 146 features on the northernmost 25% of 
the Venus surface have been classified as impact craters [Basilevsky et al., 1987]. 
These features range in diameter from approximately 10 to 140 km. A Venera 
image of the largest identified impact crater, Klenova, is shown in Figure 3.1. 

The density of impact craters observed by Venera 15/16 may be used to 
estimate the mean crater retention age of the northern quarter of the Venus surface. 
A key and highly uncertain parameter in such estimates is the crater production 
rate on Venus and its variation with time. Different assumptions regarding this 
parameter can lead to quite different values for the mean surface age. Ivanov et ai 
[1986], following Hartmann [1977], assumed that the cratering flux on Venus has 
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been equal to that for the Moon since the end of lunar mare emplacment. With 
such an assumption the mean crater retention age is about 1 b.y. Schaber et al. 
[1987a], in contrast, suggested that the cratering flux for Venus is closer to that of 
the Earth, which is nearly triple the lunar cratering flux [Grieve, 1984]. A crater 
retention age of only a few hundred million years results. 

Topographic profiles of a number of impact craters were obtained by the Venera 
15/16 radar altimeters [Ivanov et al., 1986], Basilevsky et al. [1987] give apparent 
depths (measured from rim crest to crater floor) for 20 craters. These features, with 
diameter D ranging from about 25 to 140 km and depth varying from 200 to 2400 
m, are depicted on a depth-diameter plot in Figure 3.2. 

The crater topographic data are subject to several sources of error. The first is 
the uncertainty in the measurement of the average elevation over the altimeter 
footprint. At ± 50 m this uncertainty is small compared with the relief of large impact 
craters on Venus [Barsukov et al., 1 986]. A second source of error is the smoothing 
effect of the 50-km-diameter footprint of the altimeter [Barsukov et al., 1986]; this 
smoothing yields an underestimate of topographic relief. Examination of the 
apparent depth of the 20 craters in Figure 3.2 suggests that depths are significantly 
underestimated only for craters with diameters less than or equal to about 30 km. 
We do not consider further the topography of craters less than about 30 km in 
diameter in order to avoid this source of bias. A third potential error is the mis- 
identification of features as impact craters in the tabulation of Basilevsky et al. 
[1987]. We do not include Cleopatra Patera in our analysis, because its deep inner 
floor, breached rim, and associated infilling by plains materials suggest a volcanic 
origin for this feature [Schaber et al., 1987b]. Yablochkina and Tsvetaeva are 
noticeably noncircular in map view [Basilevsky etal., 1987], and the possiblity of a 
volcanic origin for these two features has been raised [G. G. Schaber, personal 
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communication, 1987]; without further information, however, we include them in our 
analysis. After eliminating Cleopatra and smaller features with significantly under- 
sampled depths, there are 15 craters of probable impact origin with topographic 
profiles measured by Venera 15 and 16; these craters are listed in Table 3.1. 

The geographical distribution of the craters in Table 3.1 appears to be non- 
random (Figure 3.3). Two-thirds of these craters fall within a region extending over 
the pole and covering only about 10% of the Venera 15/16 survey area, or 2-3% of 
the Venus surface. A nearest-neighbor test of intercrater distances [Daw's, 1973] 
indicates that a random distribution of crater locations may be rejected with greater 
than 90% confidence in favor of a clustered distribution. This result suggests that 
there is some systematic experimental bias in the locations of impact craters with 
measured depths. These craters are only a small subset of the total observed to 
date, however, and the overall distribution or craters imaged by Venera 1 5 and 1 6 
appears to be random [Basilevsky et ai, 1 987]. These observations raise the 
question of whether the craters with known topography are globally representative. 
Available Pioneer Venus elevation data [Pettengill et a!., 1980], referenced to the 
mean radius of 6051.9 km [Ford, 1986], show that the distribution of elevations of 
the regions surrounding the craters in Table 3.1 spans the range -1 .5 to 1 .2 km and 
is given approximately by a normal distribution with mean -0.2 km and standard 
deviation 0.7 km (Figure 3.4). Further, 13 of the 15 craters occur on smooth or 
rolling plains (Figure 3.4), units widespread on the northern quarter of the Venus 
surface [Barsukov et ai, 1986; Basilevsky et ai, 1987]. We therefore assume that 
the crater relief measurements are representative of other plains units at similar 
elevations elsewhere on Venus. 
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APPROACH 

There are several possible procedures for testing viscous relaxation models 
against the preserved topography of impact craters. One approach would be to 
separate the craters with measured topography into candidate fresh and relaxed 
craters and seek viscous relaxation models that map the topography of a fresh 
crater into that of a relaxed crater of the same or similar diameter. Such a 
methodology was adopted by Hall et al. [1981] in their study of the importance of 
viscous relaxation as a crater modification process on the Moon. Unfortunately, the 
combined effects of limited horizontal resolution and the small number of craters 
with measured topography preclude such a treatment for Venus at present. 

Instead, we adopt the following approach. We employ apparent depth as the 
sole measure of crater relief on Venus, because estimates of depth are likely to be 
more robust than are the details of crater topographic profiles. From Venera data 
we calculate the mean and variance of observed crater relief on Venus as a 
function of crater diameter. From observations of lunar craters and scaling 
arguments, we then estimate the mean and variance of the apparent depth that 
fresh impact craters would display on Venus as a function of crater diameter. 
Support for this estimation procedure comes from the observation of Ivanov et al. 
[1986] that large craters on Venus are morphologically similar to impact features of 
comparable size on other planetary bodies. In particular, all of the observed impact 
craters on Venus with D > 30 km are morphologically complex, with terraced inner 
walls and central peaks on a flat floor [Basilevsky et al., 1987]. Given the two 
relations for observed and fresh crater populations we test theoretical models of 
viscous relaxation. For specific model parameters and crater diameter, we 
compare the crater depth predicted by the viscous relaxation model with the 
observed distribution of crater depths. We establish a formal statistical criterion for 
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the rejection of specific viscous relaxation models as being inconsistent with the 
observations. 

Depths of Observed and Fresh Craters on Venus 

A relationship between depth and diameter of craters on Venus may be 
obtained from the data in Table 3.1 . Linear regression of log diameter D versus log 
depth d v yields the relation 

log d v = (-0.83 ± 0.03) + (0.34 ± 0. 1 7) log D (3. 1 ) 

where both D and d v are expressed in kilometers. 

Because of the dependence of crater characteristics on the surface gravitational 
acceleration g [Melosh, 1980] and the similarity in g for Earth and Venus, terrestrial 
data are the natural choice from which to estimate initial depths of fresh craters on 
Venus. There may be systematic differences between craters on the Earth and 
Venus due to increased fallback into the crater of atmospherically decelerated 
ejecta [Settle, 1980] and increased volume of impact melt [Grieve and Head, 1982] 
on Venus, but such differences have not been detected at Venera resolution and 
appear to be minor contributions to overall relief. The relationship between 
diameter and apparent depth for complex terrestrial craters was examined by 
Grieve and Robertson [1979], who derived log d = -0.85 + 0.31 log D for craters in 
sedimentary targets. No regression was performed on the widely scattered data 
from craters in crystalline targets, but such craters were reported to be 
approximately twice as deep as structures of the same diameter in sedimentary 
targets, a result that would be consistent with the approximate relation log d = -0.55 
+ 0.31 log D. The principal source of error in estimates of apparent depth for 
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terrestrial craters is the correction for rim erosion. Examination of the available 
information on the geometry of the complex craters included by Grieve and 
Robertson [1979] suggests that initial depths are well known for no more than 
about five such structures. 

On account of the large uncertainties in the initial depths of complex terrestrial 
impact craters, we have elected to utilize the depths of 52 comparatively fresh lunar 
craters of Copernican and Eratosthenian age. The significantly larger number of 
observations will ensure that our predictions of crater depth, scaled to Venus 
conditions, will have small variance. Stereophotogrammetric and shadow- 
measurement data from Pike [1980a] and maps from the Lunar Topographic 
Orthophotomap (LTO) and Lunar Aeronautical Chart (LAC) series were used to 
obtain depths. Errors in depth measurements are tens of meters for LTOs and 
1 00-200 m for LACs. 

We assume that the initial topographic profile of a large crater on Venus differs 
from that of a fresh lunar crater of the same diameter by a simple multiplicative 
scaling factor. Observations of crater dimensions on the terrestrial planets and 
satellites have been invoked to suggest two possible scaling rules: that either the 
depth of complex craters or the crater diameter D t marking the transition from 
simple to complex crater morphology varies inversely with the surface gravitational 
acceleration g [Pike, 1977, 1980b; Grieve and Robertson, 1979; Basilevsky and 
Ivanov, 1982]. These two relations are mutually exclusive. If the depth of complex 
craters scales inversely with gravity, then d v = (g m / g v ) d m , where the subscripts m 
and v refer to the Moon and Venus, respectively. Under the alternative assumption 
that transition diameter D t scales inversely with gravity, a relation describing the 
depth of complex craters on Venus has been derived graphically from lunar data by 
Ivanov et al. [1986]. Such a relation may be given in analytic form, d v = (g m / g v ) 
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( 1 * a i) d m , under the assumptions that the ratio d/D for simple craters is independent 
of g and that the parameter ai in the relation log d ■ ao + a-i log D for complex 
craters is similar on all of the terrestrial planets. The factors for scaling the 
topography of complex craters from the Moon to Venus under the assumptions that 
d or D{ scales inversely with g are 0.183 and 0.305, respectively. Regression 
analysis of the lunar data can then be used to predict the initial depth of a complex 
crater of diameter D on Venus 

log dj = (a© ± 0.01 ) + (0.31 ± 0.04) log D (3.2) 

where ao is -0.74 under depth scaling or -0.52 under transition-diameter scaling 
and both D and dj are expressed in kilometers. The means and 95% confidence 
intervals for the mean for observed and fresh crater depths are shown in Figure 3.2; 
the fresh crater depths are those given by the depth-scaling rule. We note that the 
depth interval spanned by the two scaling rules is similar to the range of depths 
inferred by Grieve and Robertson [1979] for complex terrestrial craters in 
sedimentary and crystalline targets, although as noted above the errors in the initial 
depths of such terrestrial craters may be large. 

In addition to crater depth, we need an estimate of the complete topographic 
profile of a fresh complex crater on Venus in order to evaluate the change in 
apparent depth associated with a particular model for viscous relaxation. We used 
the profiles of three comparatively young lunar impact craters having diameters that 
approximately span the diameter range of the Venera observations: Langrenus (D 
= 132 km), Eratosthenes (D = 58 km), and Timocharis (D = 34 km). The topography 
for Langrenus was derived from Earth-based radar [S. Zisk, personal 
communication, 1986], that for Eratosthenes was obtained from LAC 58, and that 
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for Timocharis was obtained from LTO 40B3. For all three craters, an azimuthally 
averaged profile was calculated. These profiles were stretched uniformly in the 
radial direction to match specific crater diameters in the Venera data [Ivanov et ah, 
1986; Basilevsky et ai, 1987]. The elevations were adjusted in the radial stretch 
using the empirical relations of Pike [1 977] expressing the variation of key crater 
morphological parameters with diameter. For example, Pike showed that rim 
height is proportional to D 0 - 4 , so the rim height at D + AD is multiplied by (1 + 
AD/D)0* 4 i1+ 0.4 AD/D. The profiles were then multiplied by the factor 
appropriate to depth or transition-diameter scaling derived above to convert the 
lunar profiles to initial topographic profiles for Venus craters. 

Statistical Interpretation of Viscous Relaxation Models 

Given the regression curves and confidence intervals for the depths of both 
fresh and observed impact craters on Venus, we may estimate from Figure 3.2 the 
extent of viscous relaxation at each diameter such that the mean depth of partially 
relaxed craters is no longer consistent with the observations. For example, the 
predicted average depth of a fresh crater 140 km in diameter, from Figure 3.2, is 0.8 
km. The estimated lower bound at 95% confidence for the mean depth of observed 
craters on Venus at D = 140 km is 0.5 km. Therefore viscous relaxation could 
reduce the mean depth of craters at this diameter to as little as 0.5 km without 
making this datum inconsistent with the observations at the given confidence level. 
In other words, any viscous relaxation model that produced more than about a 35% 
decrease in the apparent depth of craters of this diameter for the given initial depth 
would yield craters significantly shallower than are observed. 

We may formalize this approach as follows. For a fixed D, we assume that all 
craters are of the same age, equal to half the mean crater retention age. We further 
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assume that the initial crater depths at fixed diameter are log-normally distributed, 
with mean log dj and variance of the mean Sj 2 given by linear regression on the 
scaled lunar data. Note that Sj 2 is also a function of D, with the smallest values 
near the center of the data field. Viscous relaxation will reduce the mean depth to 
8(D) (Figure 3.5), which we compare with the observed crater population, taken 
also to be log-normally distributed with mean log d v and variance of the mean s v 2 
at the same diameter. We test the null hypothesis, H 0 : d v = 8, by forming the 
Smith-Satterth waite T statistic [Devore, 1982, p. 292] 

T . (d v - J)/tf ♦ s i V ' 2 (3.3) 

with v degrees of freedom given by 

V - (8* + sf) 2 / [sj / (n y - 1) + s, 4 / (a - 1)] (3.4) 

In the limit where the variance of scaled lunar crater depths is negligible compared 
with the variance in observed crater depths on Venus, the one-sample f test may be 
used, although in practice we use the Smith-Satterthwaite test in all cases. The 
null hypothesis is rejected at the significance level a in favor of the alternative 
hypothesis, Hi : d v > 8, if T £ t^y, where t^y is the 100(1 -a) cumulative percentile 
of the t distribution with v degrees of freedom. For a value of T fixed by a specified 
amount of viscous relaxation, we can solve for the confidence level 1-a at which 
the null hypothesis may be rejected by integrating the f distribution [Selby, 1967] 



1 - a 
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where r is the gamma function. The parameter a varies from 0 to 0.5 in the 
one-taiied test; here it provides a measure of the likelihood of a particular viscous 
relaxation model given depths of observed and fresh craters. 


VISCOUS RELAXATION MODEL 

We assume that the crust and upper mantle of Venus behave as linearly 
viscous fluids. This simple rheology has successfully described the response of 
the terrestrial mantle to deglaciation [e.g M Cathles, 1975] and should be applicable 
to most of the Venus interior because of the high surface temperature and long time 
scales involved. Laboratory measurements of high-temperature creep in rocks 
typically show a nonlinear relationship between stress and strain rate [e.g., Ashby 
and Verrall, 1977], which we incorporate in an approximate fashion by means of an 
effective viscosity. Near-surface frictional sliding [Brace and Kohlstedt, 1980] is 
also parameterized by an effective viscosity. We consider below the possible 
modifying effects if portions of the crust or upper mantle of Venus display elastic 
behavior on long time scales. 

By direct numerical integration of the equations of motion of a linearly viscous 
fluid [Cathles, 1975], relaxation of topography can be evaluated for a medium with 
an arbitrary depth dependence of viscosity. Different densities in the crust and 
mantle are accommodated by solving the equations of motion in the two regions 
subject to appropriate boundary conditions. A detailed presentation of this 
treatment is given in Appendix 3.A. This method is an improvement over most 
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previous planetary applications of viscous relaxation models. In simple 
formulations involving a layer over a halfspace that include approximately the 
effects of both depth-dependent viscosity and isostatic compensation of surface 
topography [e.g., Solomon et al., 1982b], discontinuities in density (crust) and 
viscosity (lithosphere) must occur at the same subsurface horizon; this restriction is 
not required here. Further, the dependence of viscosity on thermal gradient can be 
explicitly included rather than resorting to isoviscous halfspace models [Passey 
and Shoemaker, 1981] or invoking an exponential variation with depth [Parmentier 
and Head, 1981]. 

On Earth, the differences in creep behavior between crustal and mantle 
materials can result in a weak lower crust separating stronger crustal and upper 
mantle layers. This situation arises for certain combinations of crustal thickness and 
thermal gradient, including those appropriate to many continental regions. Evi- 
dence supporting this view includes the aseismic nature of the lower continental 
crust [Chen and Molnar, 1983], the narrowness of the zone of strain accumulation 
across major strike-slip faults [Turcotte et al., 1984], and the two scales of defor- 
mation in zones of continental extension [Zuber et al., 1986]. The crust-mantle 
boundary, therefore, is a major mechanical discontinuity within the lithosphere. By 
analogy with terrestrial behavior, we consider a model for Venus consisting of a 
chemically and mineralogically uniform crust overlying an ultramafic mantle to 
describe the principal compositionally-induced variation in rheology within the 
upper few tens of kilometers. 

The effective viscosity in the crust and mantle depends strongly on temperature. 
High-temperature deformation of planetary materials, in general, can be described 
by a flow law of the form 
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e = A o n exp(-Q/RT) (3.6) 

where e is a strain rate, a is a corresponding deviatoric stress, R is the gas 
constant, T is the absolute temperature, and A, n, and Q are empirically derived 
constants. For linearly viscous flow under uniaxial compression (the state of stress 
under which creep measurements are typically made), the consititutive relation is a 
= 2 r| e, where rj is the dynamic viscosity. Substituting from (6), the effective 
viscosity has the form 


t](T) = [oO- n >/2A] exp (Q/RT) (3.7) 

For the mantle, we use the flow law for polycrystalline olivine (Fo 9 i) determined by 
Goetze [1978]. For the crustal layer we use the flow law for Frederick diabase 
given by Caristan [1982]. Shelton [1981] also reported a flow law for this material; 
the predicted strain rate is about an order of magnitude lower than that given by the 
relation of Caristan [1982] over the range of temperatures and stress differences 
relevant to the viscous relaxation problem. The difference may be the result of 
different experimental procedures (particularly the nature of the confining medium) 
or a greater proportion of pyroxene and a lesser proportion of weaker feldspar in 
the particular samples used by Shelton. The composition of the Venus surface 
measured at the Venera 13 and 14 landing sites [Surkov et al., 1983] is closer to 
that of the diabase samples of Caristan than that of Shelton's samples, so 
Caristan's flow law is preferable for our application. However, viscous relaxation 
calculations were also carried out using Shelton's results as a test of sensitivity to 
variations in crustal rheology. 

Creep rates in rock can be strongly affected by the presence of water, 
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particularly In the form of thin films between silicate grains fe.g. , Ashby and Verrall, 
1977). Both Caristan [1982] and Shelton [1981] endeavored to dry their diabase 
samples prior to measurement of strain rate, and the olivine analyzed by Goetze 
[1978] contained only a small (0.3 %) weight fraction of water. Nonetheless, 
neither the relative water contents of the laboratory samples and the Venus interior 
nor whether any differences significantly affect the flow laws are known. 

We assume that the temperature increases linearly with depth from an average 
surface value, equal to the surface temperature at the mean planetary elevation 
[Seiff et al., 1980]. Because the dependence on stress observed in laboratory 
samples is nonlinear, the effective viscosity must be evaluated at some character- 
istic stress difference o, which we take to be the average initial stress difference 
within a cylindrical volume of radius D and depth D/2 centered beneath the crater. 
These dimensions are somewhat arbitrary but approximately reflect the volume 
within which significant deviatoric stress occurs. In practice, the characteristic 
stress is solved iteratively by specifying an initial value for c, calculating the 
average initial viscous stress deviator within the designated volume, substituting 
this new value for o, and repeating the calculation until convergence occurs. 
Because such calculations would be extremely time consuming if conducted for 
every model, we determine 0 in this fashion at each crater diameter only for a 
single reference model with representative parameter values. We discuss this 
approximation further below. 

Near the Venus surface rock strength is likely limited by frictional sliding on 
faults rather than by ductile flow. This effect can be parameterized by specifying a 
brittle-ductile transition depth above which the effective viscosity increases linearly 
with depth. The viscosity in both the brittle and ductile regions is then proportional 
to lithospheric strength as defined by Brace and Kohlstedt [1980]. Values of the 
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characteristic stress a and the depths of the brittle-ductile transitions in the crust 
and mantle are determined simultaneously for the reference model. Adopted 
model parameters are given in Table 3.2. 

Largely for numerical convenience, we impose a lower bound of 10 20 Pa s on 
the effective viscosity in either the crust or the mantle. Substantially lower values 
are not likely to be characteristic of long-term bulk deformation of the Venus 
lithosphere. The results reported below are not sensitive to the value selected for 
this lower limit. 

As an illustrative comparison, we have applied the viscous relaxation model to 
a large impact crater on the stable interior of a terrestrial continent. On Earth, 
erosion modifies crater topography at rates of order 10 m/m.y. [Grieve and 
Robertson, 1979], and a primary issue is whether viscous relaxation of relief occurs 
at comparable or greater rates and can therefore be resolved. We take D = 140 km 
as the largest preserved terrestrial impact structure [Grieve and Robertson, 1 979 ], 
and we assume that the initial topography may be obtained from that of the lunar 
crater Langrenus under the scaling that crater depth varies as g* 1 . In this test case, 
we select parameters so as to maximize the rate of relaxation. We adopt H = 35 km 
and dT/dz = 20 K/km, and we assume that the flow law for the continental crust is 
that for quartz [Brace and Kohlstedt, 1980]. We also decrease the volume over 
which the characteristic stress o is obtained to that contained within the cylinder of 
radius D/2 and depth D/4, giving a = 15 MPa and a depth of 15 km for the brittle- 
ductile transition in the crust. We find rates of relaxation of topography of a few 
meters per m.y., which result in relaxation of 45% of the crater relief in 200 m.y. 
These rates are an order of magnitude smaller than the rate of removal of relief by 
erosion. We conclude that viscous relaxation has not been an important 
modification process for impact craters on stable terrestrial continents. 
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APPLICATION TO VENUS 

The key parameters in the viscous relaxation model that are of geological and 
geophysical interest are the crustal thickness H and the average lithospheric 
thermal gradient dT/dz (see Appendix 3.A). Comparison of the observed relief of 
impact craters with that predicted from the model can be used to derive upper 
bounds on these parameters. Since viscosity decreases with increasing temper- 
ature, relaxation will occur more rapidly at higher temperature, or at a higher value 
of dT/dz for a fixed surface temperature. Because the crust is less viscous than the 
mantle at a comparable temperature, topographic decay will be more rapid for a 
thicker crust. We can therefore anticipate that a trade-off exists between dT/dz and 
H when other parameters are fixed: a given amount of topographic relaxation in a 
specified time can be accommodated by a low thermal gradient and a thick crust or 
alternatively by a higher thermal gradient and a thinner crust. In terms of the 
statistical approach discussed in the last section, models with higher thermal 
gradients and thicker crusts can be rejected more easily than models with lower 
thermal gradients and thinner crusts. We adopt the a = 0.05 significance (95% 
confidence) level for the statistical test of upper limits to dT/dz and H. The risk of 
incorrectly rejecting a (dT/dz, H) combination (type I error) is therefore only 5%. The 
corresponding probability of incorrectly accepting a (dT/dz, H) pair (type II error) is 
about 20% [ Devore , 1982, p. 246]. Since we seek upper limits to dT/dz and H, the 
possibility of a type I error is more serious and should be minimized in spite of the 
higher probability of the alternative type of error. 

To ensure that the derived upper bounds on dT/dz and H are valid, we should 
utilize observations and adopt modeling approaches such that errors are likely to 
lead to overestimates of these parameters. For instance, observed crater depths 
are likely to be underestimates of the true depths because of the smoothing effect 



52 


of the altimeter footprint (see above), and because such additional processes that 
degrade crater relief as volcanism or erosion have been neglected. Under these 
conditions the amount of viscous relaxation will be overestimated, which will raise 
the upper bounds on dT/dz and H. The possibility that young impact craters on 
Venus are systematically shallower than terrestrial or scaled lunar craters because 
of enhanced ejecta fallback [Settle, 1980], increased impact melt production 
[Grieve and Head, 1982], or more rapid initial strain rates during the transient 
decay of impact heat [O’Keefe and Ahrens, 1977] from a target at higher ambient 
temperatures, will also yield overestimates of the extent of any long-term viscous 
relaxation. The characteristic stress, described in the previous section, is also 
deliberately underestimated in order to overestimate the effective viscosity for a 
specified material and temperature; this in turn yields larger upper bounds on dT/dz 
and H. 

Beyond crustal thickness and average thermal gradient the remaining para- 
meters in the viscous relaxation problem are the crater age t, the initial state of 
isostatic compensation c (see Appendix 3.A), and the choice of a scaling rule to 
estimate the depths of fresh impact craters from the depths of lunar craters. We will 
consider average crater ages of 50 m.y. and 500 m.y., corresponding to mean 
crater retention ages of 100 m.y. and 1 b.y., respectively; these values approx- 
imately span the range in reported values for the northern 25% of the Venus 
surface [Ivanov et al., 1986; Schaber et al., 1987a; Grimm and Solomon, 1987]. 
The parameter c is taken to be 0 or 1 , corresponding to zero or complete initial 
local compensation of crater topography. For all calculations described below we 
shall assume that crater depth scales as g _1 . The alternative scaling law in which 
transition diameter scales as g* 1 yields greater depths for fresh craters at each 
diameter. This greater relief results in a greater effective stress beneath the crater, 
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which in turn drives more rapid viscous relaxation. However, the accumulated 
fractional change in relief before the depth of a partially relaxed crater becomes 
inconsistent with observed depths is also larger for this scaling law (about 60% at 
140 km diameter instead of 35%). These two effects tend to offset one another, so 
that the final distribution of confidence levels is not very different from those derived 
under the adopted scaling rule. 

Initially Uncompensated Craters 

The viscous relaxation of a crater 140 km in diameter is depicted in Figure 3.6. 
This diameter corresponds to the largest crater on Venus with a measured depth 
(Table 3.1 and Figure 3.2). Relaxation is shown for two sets of values of thermal 
gradient and crustal thickness; for both sets of calculations the crater is assumed to 
be initially uncompensated. With a low thermal gradient (10 K/km) and a thin crust 
(10 km) there is relatively little relaxation of relief even in 500 m.y. (Figure 3.6a), 
whereas with a somewhat higher thermal gradient (15 K/km) and a somewhat 
thicker crust (15 km) relaxation is nearly complete after only a few tens of millions of 
years (Figure 3.6b). 

A similar set of calculations for a crater 32 km in diameter is shown in Figure 
3.7. This diameter corresponds to the smallest Venus crater for which the mea- 
surement of depth is likely to be reliable (Table 3.1 and Figure 3.2). Compared 
with the case for D = 140 km, there is a comparable relaxation of relief for a thin 
crust and low thermal gradient but substantially less relaxation for a thicker crust 
and higher gradient. 

This behavior may be understood by consideration of the constant x governing 
the exponential decay with time of surface topographic relief; these time constants 
are shown in Figure 3.8 for the four models in Figures 3.6 and 3.7. (Formally, there 
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are two time constants t-| and %z for each model, but the term corresponding to the 
larger time constant x-j makes only a small contribution to the topography in this 
case; see equation A40 in Appendix 3. A). For given dT/dz and H, values of x at 
different D differ by a multiplicative constant; the lower stress difference beneath a 
smaller crater results in a higher effective viscosity, which in turn leads to a longer 
relaxation time. The largest contributions to the Hankel transform of the crater 
topographic profile are at wavelengths between 0.5 and 3 D. For models with 
dT/dz « 15 K/km and H = 15 km, the time constants at these wavelengths for a 
crater with D = 32 km are of order 10 8 yr, whereas those for a crater with D = 140 
km are of order 1 0 5 - 1 0 7 yr. Much of the relief of the smaller crater therefore 
remains after 500 m.y. (Figure 3.7b) while the larger crater relaxes completely over 
the same time (Figure 3.6b). 

The situation differs for models with dT/dz = 10 K/km and H = 10 km. The time 
constants at the relevant wavelengths for the crater with D = 32 km are near a focal 
minimum in x controlled by the low viscosity zone in the lower crust. The important 
time constants for the crater 140 km in diameter, however, are near a local maxi- 
mum in x controlled by the high viscosity upper mantle. Because of the stress- 
dependence of effective viscosity, the values of x near the respective extrema are 
similar and of order 10 9 yr. At a given time after crater formation, the fractional 
relaxation for craters of both diameters is similar under these circumstances 
(Figures 3.6a and 3.7a). 

We have conducted calculations of the sort shown in Figures 3.6 and 3.7 for a 
wide range of values of both dT/dz and H. For a crater 140 km in diameter, the 
fractional relief remaining after 50 m.y. is shown as a function of these two para- 
meters in Figure 3.9a. The trade-off between dT/dz and H is apparent. Corres- 
ponding confidence intervals for rejecting parameter values are given in Figure 
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3.9b. Figure 3.10 shows the remaining relief and the associated confidence levels 
for parameter rejection for a crater 140 km in diameter and 500 m.y. in age. For 
given values of dT/dz and H, the confidence level for model rejection is higher at a 
500 m.y. age than a 50 m.y. age, because the larger predicted amount of relaxation 
of topographic relief is more difficult to reconcile with the observed crater depths. 
The results for a similar suite of calculations for a crater 32 km in diameter are 
shown in Figures 3.11 and 3.12. At this crater diameter the fractional relief 
remaining is insensitive to crustal thickness for H > 20 km because the upper 
mantle does not play a significant role in the flow accommodating relaxation of 
shorter wavelengths. 

A summary of 95% confidence limits to paired values of dT/dz and H as 
constrained by these models with initially uncompensated crater relief is given in 
Figure 3.13. For D = 140 km, lithospheric thermal gradients in excess of 20 K/km 
are possible only when the crust is thinner than 10 km and then only if most craters 
are young (of order 50 m.y.). If the average lithospheric thermal gradient is as low 
as 10 K/km, then a 15-20 km thick crust is possible, again only if most craters are 
young. However, if the age of the Venus surface is comparable to the largest 
current estimates (mean crater ages of order 500 m.y.), then acceptable models are 
restricted to H £ 10 km and dT/dz £ 10 K/km. The range of possible combinations 
of H and dT/dz is larger for D = 32 km, since relaxation is generally slower for 
shorter wavelengths. For such features, the crustal thickness is unconstrained 
provided that the observed craters at this diameter are young and that dT/dz < 20 
K/km. For a 500 m.y. mean crater age, similar behavior independent of H requires 
dT/dz <10 K/km. 

We have also examined viscous relaxation for a Venus crater with D = 50 km, a 
value near the median diameter of craters in Table 3.1. Although the relaxation of 
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a 50-km-diameter crater is somewhat closer to that for D = 32 km than that for D = 
140 km, the inferred limits to crustal thickness and thermal gradient are closer to 
those imposed by the relief of the largest craters, with upper limits to H of 20, 10, 
and 5 km for dT/dz = 10, 20, and 30 K/km respectively (Figure 3.13). This result is 
because the constraint on the mean depth of observed craters is strongest in the 
center of the diameter distribution than at either extreme (Figure 3.2) and so less 
viscous relaxation is allowed at a specified confidence interval. 

For a diabase crust that deforms according to the flow law derived by Shelton 
[1981], limits on dT/dz and H are less restrictive than those using Caristan's [1982] 
flow law, probably because of the greater proportion of stronger pyroxene in 
Shelton’s samples. In general, the range in allowable H is 5-10 km greater at fixed 
dT/dz than shown here. 

Initially Compensated Craters 

We have repeated these calcuations under different assumptions concerning 
the initial state of isostatic compensation of crater topography. Initial compensation 
will lessen the relaxation rate, principally because of lesser stress differences to 
drive flow. Partial to complete initial compensation of relief, therefore, is permissive 
of a broader range of acceptable combinations of H and dT/dz. Curves showing 
95% confidence limits for rejecting combinations of values for dT/dz and H for 
initially compensated craters are shown in Figure 3.13. If the characteristic stress 
difference a is taken to be equal to that for the initially uncompensated crater of the 
same diameter (Table 3.2), then the solutions at dT/dz < 20 K/km are similar for 
compensated and uncompensated cases because relaxation is limited by the large 
viscosity of the upper mantle. At dT/dz > 20 K/km, H is still restricted to be less than 
about 1 0 km. 



57 


Calculation of the characteristic stress difference 6 is much more sensitive to 
variation of dT/dz and H for initially compensated relief (c = 1) than for the situation 
with c = 0. For an initially uncompensated crater 1 40 km in diameter, a varies over 
the narrow interval 2.5 - 5 MPa for the ranges of dT/dz and H shown in Figure 3.13, 
whereas for c = 1 , o varies over the interval 0.1 - 3.5 MPa. For a reference model 
given by dT/dz = 13 K/km and H = 1 3 km, o is 4 MPa and 0.5 MPa for c = 0 and c = 
1 , respectively. The latter value for the characteristic stress difference is probably 
artificially low, because of the large averaging volume, yielding a high effective 
viscosity and conservative upper limits to H and dT/dz for initially compensated 
large craters. For this situation, thermal gradients in excess of 20 K/km and a thick 
crust (H > 30 km) are permissible parameters if the crater is young (50 m.y.). For an 
older crater H must be less than about 20 km at dT/dz = 25 K/km. The difference in 
o between initially compensated and uncompensated cases is minor at 32 km 
diameter, because the surface and crust-mantle boundary are effectively 
decoupled at the relevant wavenumbers for H > 10 km. 

Finite Lithospheric Strength 

We have also considered the support of crater relief by an elastic-plastic 
lithosphere; see Appendix 3.B. Because the strength in the lower portion of the 
lithosphere is limited by temperature-dependent ductile flow, a lower bound on the 
extent of flexural support of relief may be converted to an upper bound on mean 
lithospheric thermal gradient. For craters 140 km in diameter the limits on dT/dz if 
relief is supported flexurally (Table 3.B1) are similar to those from viscous 
relaxation models with no initial isostatic compensation and with comparatively 
young crater ages (of order 50 m.y.) for the cases of both thick and thin crust. The 
topography of smaller craters, however, can be supported elastically in a 
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lithosphere having a crust of arbitrary thickness for dT/dz up to 30 K/km at D = 32 
km and dT/dz up to 1 8 K/km at D = 50 km (Table 3.B1 ). If an elastic-plastic model is 
appropriate to the lithosphere of Venus, the depths of smaller craters do not 
provide useful information on crustal thickness or thermal gradient. 

DISCUSSION 

The results of the viscous relaxation models presented in the previous section, 
and summarized in Figure 3.13, provide limiting values for crustal thickness and 
lithospheric thermal gradient consistent with the preservation of topographic relief 
in large impact craters on Venus. In this section we compare these limiting values 
to estimates derived independently by others, and we explore the implications of 
these results for the evolution of Venus. 

As noted above any finite strength of the Venus lithosphere can contribute 
significantly to the support of topographic relief for craters with D ^ 50 km. The 
largest (D > 100 km) craters, therefore, provide the most stringent and robust 
constraints on crustal thickness and lithospheric thermal gradient. The small 
number of such craters with measured relief makes uncertain any extrapolation to 
regional or global conditions. It should be recognized, however, that the limits on 
the extent of relaxation are derived from the entire data set depicted in Figure 3.2; 
i.e., the variance in the mean depth at large diameter is a function of the variance in 
the data at all diameters. Because of the possibility of partial elastic support for 
even the largest craters, and to maintain a conservative approach, we shall use the 
limiting values obtained under the assumption of a young mean surface age 
(Figure 3.13a) in the following discussion. 

According to Figure 3.13a, the crustal thickness beneath the largest (D > 100 
km) craters with measured depths on Venus and, by inference, beneath geo- 
logically similar plains units at comparable elevations, is less than about 10 km for 



59 


a mean lithospheric thermal gradient of 20 K/km and less than about 20 km for a 
mean lithospheric thermal gradient of 10 K/km. These values may be compared 
with other estimates of lithospheric thermal gradient and crustal thickness on 
Venus. Assuming that Venus loses heat at the same rate per unit mass as the 
Earth, the global rate of heat loss is 3.4 x 10 13 W and the average heat flux is 74 
mW/m 2 [Solomon and Head, 1 982]. Alternatively, if radioactive heat sources are 
present in Venus at chondritic abundances, then the average heat flux may be as 
low as 50 mW/m 2 [Phillips and Matin, 1984]. If all of this heat flow is transported 
through the lithosphere by conduction and if the thermal conductivity of the 
lithosphere is similar to that on Earth (3.1 W/m-K, [Parsons and Sclater, 1977]), 
then the mean lithospheric thermal gradient is 16-24 K/km. Locally the thermal 
gradient may be as low as that in terrrestrial continental shields and old ocean 
basins, about 10 K/km [Sclater et ai, 1980]; in areas of active rifting and volcanism, 
of course, the local gradient may be considerably greater than 25 K/km. 

Upper bounds to the crustal thickness on Venus of one hundred to several 
hundred kilometers have been suggested on the basis of the projected depth of the 
basalt-eclogite stability field [Anderson, 1 980] or the apparent depth of isostatic 
compensation of long-wavelength topography [Phillips et ai, 1981; Phillips and 
Malin, 1984]. These bounds may greatly exceed the actual crustal thickness if the 
crust does not fill the volume of the interior lying within the basalt stability field or if 
mantle dynamical effects contribute significantly to the support of long-wavelength 
topography. Much lower values of maximum crustal thickness were derived by 
Zuber [1987] and by Banerdt and Golombek [1988] from the observation that there 
are often two characteristic wavelengths of tectonic features in extensional and 
compressive terrains on Venus. These workers interpret this observation as 
evidence for a layered lithosphere: strong upper crustal and mantle layers 
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separated by a weak lower crust. These layers arise from the different 
temperatures at which ductile flow limits strength in crustal and mantle material 
(see Appendix 3.B for a related discussion); the observed scales of deformation 
can be quantitatively related to crustal thickness and lithospheric thermal gradient. 
Zuber [1987] found that the crustal thickness is less than about 15 km if dT/dz = 25 
K/km and less than about 30 km if dT/dz as 10 K/km. These values are quite similar 
to those obtained in this study (Figure 3.13a). 

There are two scenarios that would invalidate the constaints on crustal thick- 
ness and thermal gradient implied by Figure 3.13, but both of these possibilities are 
unlikely. First, the largest craters on Venus with measured relief could be very 
young, i.e., considerably younger than 50 m.y., so that viscous relaxation would be 
minimal for any plausible model. The ages of 40-2000 m.y. for the largest (D - 100 
to 140 km) terrestrial craters [Grieve and Robertson, 1979] do not support this view, 
however. Further, the distribution of crater morphological types, or states of preser- 
vation, appears to be independent of crater diameter [Basilevsky et al., 1987]. In 
particular, the largest craters with measured relief have intermediate states of 
preservation (Table 3.1) and therefore are inferred not to be substantially younger 
than the overall crater population. In a second scenario, the topographic relief of 
large impact craters could be completely compensated by variations in crustal 
thickness; as noted above, a crustal thickness of 30 km or more would then be 
compatible with the observed relief and inferred surface age even in the presence 
of a lithospheric gradient in excess of 20 K/km. Variations in the Bouguer gravity 
anomaly across large (D ~ 100 km) terrestrial impact structures in continental crust, 
however, are modest (~10 mgal) and do not indicate significant Moho relief 
[Masaitas et al., 1976; Sweeney, 1978]. It is doubtful that in crust of comparable 
thickness on Venus the craters in this size range should be completely 
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compensated. 

We conclude that the thickness of the crust beneath the largest (D > 100 km) 
craters, and possibly beneath plains units over at least the northern quarter of the 
Venus surface, is no more than 10-20 km thick. If this result holds for plains units at 
similar elevations elsewhere on the planet, we may obtain a simple estimate of the 
total crustal volume from the hypsometry of the Venus surface by assuming that 
long-wavelength topographic variations are locally compensated by crustal thick- 
ness variations. Measured gravity anomalies indicate that other compensation 
mechanisms also operate [Phillips et al., 1979; Phillips and Malin, 1984], so that 
the crustal volume estimated under this assumption is an upper bound. If H 0 is the 
crustal thickness beneath plains units of elevation z 0 , then the crustal volume is 
given by 


2 max 

V - J [H 0 + (2 - zj (1 + p c /Ap)] dA (z) (3.8) 

2 nin 

where dA(z) is the area with elevation between z and z + dz, and where z m j n and 
z max are the are the minimum and maximum elevation. If the Pioneer Venus 
hypsometry curve [Masursky et al., 1980; Ford , 1986] is representative of the 
planet, and taking p c = 2.9 Mg/m 3 , Ap = 0.5 Mg/m 3 , and z 0 = 6051.5 km (the 
elevation of the crater Cochran), we calculate V = 0.5 - 1 .2x1 0 10 km 3 for H 0 equal 
to 1 0-20 km. Such a volume is equivalent to a global mean crustal thickness of 
10-25 km, and is comparable to the Earth’s present crustal volume. 

Whether this volume of crust represents the total volume of crustal material 
produced on Venus is uncertain. Assuming that the present terrestrial rate of plate 
creation [Parsons, 1982] is approximately representative of the past 4 b.y., the total 
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volume of crust generated on Earth has been ~ 10 11 km 3 . If Venus has produced a 
comparable crustal volume, then recycling of older crust, by subduction, founder- 
ing or remelting, must have occurred. If 10 10 km 3 is the full time-integrated volume 
of crustal material on Venus, then the average rate of crustal generation over the 
last 4 b.y would be 2 km 3 /yr. This value coincides with the upper bound on the 
global rate of volcanic resurfacing obtained from the density of preserved impact 
craters [Grimm and Solomon , 1987]. It should be noted that the latter figure is 
insensitive to the contribution to crustal generation by igneous intrusions. Since 
the rate of intrusion of new crustal volume may be several times the rate of 
extrusion of new material, a value of 10 10 km 3 would for the crustal volume would 
suggest that the average volcanic flux on Venus in the absence of crustal recycling 
could be considerably less than 2 km 3 /y. 

CONCLUSION 

The depths of large (D > 100 km) impact craters on Venus limit the extent of 
viscous relaxation of topographic relief that has occurred since crater formation. 
Under the assumption that the crust-mantle boundary is the principal mechanical 
discontinuity of compositional origin within the lithosphere, we have developed a 
new solution to the viscous relaxation problem that incorporates a density contrast 
at that boundary and an arbitrary vertical variation in viscosity. By means of a 
statistical formalism for testing relaxation models against observed crater depths, 
we have derived trade-off curves for acceptable combinations of crustal thickness 
and lithospheric thermal gradient beneath such craters. The preservation of 
significant crater topographic relief may be achieved by a low thermal gradient or a 
thin crust. Because average thermal gradients are not likely to be less than 10 
K/km on a regional basis, a thin (1 0-20 km) crust is the most plausible explanation 
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for the degree of preservation of relief of large impact craters on Venus. This is a 
conservative estimate, as we have consistently assigned limits to mean crater age, 
apparent crater depth, amount of viscous relaxation, effective stress, and finite 
lithospheric strength in such a manner as to yield upper bounds on local crustal 
thickness and lithospheric thermal gradient. 

The large craters on Venus with measured depths, although few in number, are 
representative of a larger number of imaged craters in terms of regional elevation, 
geological unit, and state of preservation of morphological detail. It is therefore not 
unreasonable to infer that the crust does not exceed 10-20 km beneath much of the 
lowlands and rolling plains throughout the northern quarter of the Venus surface. 
Global extrapolation of this result yields a total volume of crust on Venus of about 
10 10 km 3 , comparable to the present crustal volume on Earth. This similarity, 
despite the role of plate tectonics in recycling terrestrial crust, suggests either that 
the rate of crustal generation on Venus over the past 4 b.y. has been an order of 
magnitude less than on Earth or that some form of crustal recycling has occurred 
on Venus. 

The conclusions reached in this work are sensitive principally to the accuracy of 
experimentally derived flow laws and their applicability to crustal and mantle 
material on Venus and to the small number of measurements to date of the 
topography of Venus impact craters. With further measurements of rock defor- 
mation and the improved coverage and resolution of altimetry and imaging of the 
surface of Venus to be obtained by the Magellan mision, the approach followed in 
this paper should lead to better constraints on the structure of the crust and upper 
mantle of Venus. 
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APPENDIX 3.A: 

Numerical Treatment of the Viscous Relaxation of Topography 

In this appendix we derive a numerical formulation for the gravity-driven relax- 
ation of a topographic load at the surface of a density-stratified fluid having an 
arbitrary vertical variation of viscosity. In principle a fluid having any number of 
constant-density layers can be treated by this method; the present analysis is 
confined to a two-layer model with distinct crust and mantle. 

The model is illustrated in Figure 3.A1. The vertical coordinate z is positive 
downward. The crustal thickness is H, the crustal density is p c , the mantle density 
is p c + Ap, and the viscosity is ri(z). 

Let the surface topography h s be azimuthally symmetric, that is, a function of 
horizontal distance r and time t only. In the wavenumber domain the topography 
can be represented by the zero-order Hankel transform 

oo 

C<N>-jh i (r,t)J 0 <kr)fdr (3.A1) 

0 

where J 0 is the zero-order Bessel function and k is the wavenumber. The 
corresponding inverse transform is 


h s - J h i 0) M J o < kr > k dk (3.A2) 

0 

A similar relation exists for relief h m (r,t) on the crust-mantle interface. The initial 
relief at z = H can be obtained from the initial surface relief by an assumption 
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regarding the initial state of local isostatic compensation 



(k,0) = 


CPe 

Ap 



(k,0) 


where c is the initial degree of isostatic compensation (c = 0 corresponds to uncom- 
pensated topography; c = 1 corresponds to complete Airy isostasy). 

The rest of this discussion treats the relaxation of a single harmonic component 
of relief h s at the surface and h m on the crust-mantle boundary, and the wave- 
number dependence is implicit. The complete solution for the decay of arbitrary 
surface topography is obtained by integrating over all wavenumbers via (3.A2). 

Consider the interior of an incompressible fluid away from any density discon- 
tinuities, and assume that inertial stresses are much smaller than viscous stresses. 
Then the equation of static equilibrium 


v« a + t = o 


(3.A3) 


applies, where G is the stress tensor and f is the body force per unit volume. 
Normal stress is taken to be positive in extension. The continuity condition is 

V* v = 0 (3.A4) 


where v is the fluid velocity. With the assumption of circular symmetry, expressions 
(3.A3) and (3.A4) in cylindrical polar coordinates are [Fung, 1977] 
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^ o r +|( W +a *®«“ 0 (aA5) 

9 z <5 zz + 3 r <5 f 2 + J a a + P9 - 0 (3.A6) 

3 r v, + lv + a 2 v z = 0 (3.A7) 

where we have assumed that f(r, 0, z) = (0, 0, pg) and where 9 r and 9 Z denote the 
operators 9/9 r and 9/9z, respectively. The constitutive relation for an 
incompressible viscous fluid is 


a = -p I + 2 ti D (3.A8) 

where I is the identity matrix, p is the total pressure, and D is the strain-rate tensor. 
Under azimuthal symmetry (3.A8) becomes 


= -p + 2 t| 9 f v f 

(3.A9) 

v. 


= -p + 2 r] - 

(3.A10) 

= -P + 2 ti 9 z v z 

(3. All) 


°rz= T >< 3 , v z + 3 z v ,> 


(3.A12) 
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re 



(3.A13) 


Boundary conditions require that c t 2 and o zz vanish at the free surface; that v r> 
v z , c a , and a zz are all continuous across the crust-mantle boundary; and that v r and 
v z vanish at large depth. Since the radial variation of these quantities will be given 
by a superposition of Bessel functions, we seek a system of equations which 
describes their vertical variation. We therefore must transform the equilibrium 
equations and constitutive relation so that the only derivatives appearing are with 
respect to z. We make use of the following properties of Hankel transforms 
(denoting zero- and first-order transforms by the superscripts (°) and CO, 
respectively): 


[d r f(k)] (0) 




(3.A14) 


[3 r W] ( ” 



(3.A15) 


| (V + 7 ■ T f ) J i < kr > r * - - k2 ,m (3.A16) 

n t 


where J-j is the Bessel function of order one. Transformed versions of equations 
(3.A12) and (3.A7) can be written, respectively, as 


^ J 1 ) (i) M 

3 z V r ■®B / n 


+ kv 


( 0 ) 


(3A17) 
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a 2 v[ 0) = -kv r 0) (3.A18) 

Using (3.A7) to solve (3.A1 1) for p, and substituting into (3.A9) and (3.A10), we can 
express the remaining nonzero stress components in terms of v r and a zz : 

v. 

CT rr * a 22 + 4t 1 a r v r + 2r \ J ( 3 ' A1 9 ) 

V. 

a ee = °zz + 2l l a r v r + 4l l 7 (3.A20) 

Substituting (3.A19) and (3.A20) into (3.A5) gives 

3 , + 4 n ( 3 „ v , + f 3 , v, - t) = ° (3.A21 ) 

r 

which, from (3.A15-3.A16), becomes 

d z o%> = k + 4 ii k 2 vj ,) (3.A22) 

in the wavenumber domain. 

The incorporation of an equation for d z ct zz requires a modification of the 
definition of stress. Since the hydrostatic component of pressure makes no 
contribution to the flow, we can express the viscous stresses relative to this state. 
Write the total pressure as 


P = Po + Pi 


(3.A23) 
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where p 0 = pgz and p-j is the perturbation due to flow. Then (3.A1 1 ) becomes 

O77 = O77 ■ pgz (3.A24) 

where o^’ = - p 1 + 2q 3 Z v z is the vertical normal stress relative to the hydrostatic 
pressure. The shear stress is unaffected, so a a ' = g u , and (3.A6) becomes 

d z°'zz +d r C 'rz + J C rz = 0 ( 3A25 ) 

For convenience, we drop the primes, and measure all stresses relative to the 
hydrostatic state unless otherwise noted. In the wavenumber domain, (3.A25) is 

<£' = - k <£’ (3.A26) 

We may combine (3.A17), (3.A18), (3.A22), and (3.A26) into the system of 
equations 


[v (1) l 






r v (i) i 

V r 


0 

k 

1/T| 

0 


V r 

(0) 

V z 


-k 

0 

0 

0 


(0) 

v z 

rz 


4r|k 2 

0 

0 

k 


o (1) 

rz 

o {0) 

ZZ 


0 

0 

-k 

0 


o (0> 

ZZ_ 


(3.A27) 


We will call to = [v/ 1 ), v z (°), 0^0), c zz (°)] T the velocity-stress vector, and we drop the 
Hankel transform superscripts for convenience. An analog to equation (3.A27) 
holds in two-dimensional cartesian coordinates with co = [iv x , v z , ia xz , o zz ] T , where 
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all quantities are Fourier transformed [Cathles, 1975]. 

Equation (3.A27) satisfies continuity of velocity and stress everywhere in the 
fluid away from density discontinuities. At such boundaries, continuity of velocity 
and stress must be ensured on the deformed interface. However, if the boundary 
relief h is small compared with the wavelength and the layer thickness, then to first 
order the boundary conditions may be applied on the undeformed interface, 
requiring only a correction Ap g h to a zz [Solomon et al., 1982a; Richards and 
Hager, 1984]. 

The boundary conditions at the surface are 


<%(0) - 0 

(3.A28) 

°zz(°) - Pc 9 h s 

(3.A29) 


The boundary conditions at the crust-mantle interface are implicitly satisfied by 
(3.A27), with the exception of the correction to normal stress 

Aob(H) = -Apgh m (3.A30) 

The final two boundary conditions, that the velocity components vanish at large 
depth, are also implicitly satisfied by (3.A27). 

Evaluation of (3.A27-30) is more convenient if these equations are expressed in 
nondimensional form, which is achieved by dividing each dimensional quantity by 
an appropriate reference value. We choose H as the reference length, p c as the 
reference density, and rj 0 =tj(0) as the reference viscosity. It follows that the 
reference values for time, velocity, and stress are t 0 = q 0 /(p c g H), v 0 = p c g H 2 /rj 0 , 
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and c 0 = p c g H, respectively. Equations (3.A27) and (3.A28) are unchanged with 
nondimensional quantities, whereas (3.A29) and 3.A30) become o 22 (0) = h s and 
Ao zz (H) = -Ap h m , respectively. 

A linear combination of three independent solutions to (3.A27) is required to 
satisfy the three conditions (3.A28-A30). Let these three solutions, in non- 
dimensional form, be ©('), i = 1, 2, 3. Two of these solutions, coCO and to( 2 ), are 
specified at depth do » H, and (3.A27) is numerically integrated to z = 0 using a 
fourth-order Runge-Kutta method [Shampine et ai, 1976]. Although any two 
independent vectors are sufficient, the most accurate results ensue when the 
specified values of these vectors have the form [A, B, A, B] T , where A and B are 
constants [Cathles, 1975]. Accordingly, we set co< 1 ) (do) = [1,0, 1 , 0] T and co( 2 )(do) = 
[0, 1 , 0, 1] T . The effect of the crust-mantle discontinuity is ignored in the integration 
of (3.A27) with these vectors, so the third solution must represent the apparent 
normal stress discontinuity at z = H. This requires that co( 3 ) (H) = [0, 0, 0, 1] T . Again, 
the upper limit of integration is z = 0. 

We now form linear combinations of these solutions to match the boundary 
conditions. Denoting ©jj as the j 1h component of the nondimensional velocity-stress 
vector g)(0, the boundary conditions (3.A28)-(3.A30) become 


oJO) = I a. ©.3(0) = 0 

(3.A31) 

a Z2 (°) = | a. ajj 4 (0) = h s 
i«1 

(3.A32) 

Ao^fH) = <* 3 io 34 (H) = -Ap h m 

(3.A33) 


Equation (3.A33) allows immediate solution for a 3 = -Ap h m , since co 34 (H) = 1 . 
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We are left with a 2 x 2 system for ai and a 2 : 


© 13 (°) 

© 23 (°) 


«i 


A P h m“te(°) 

®, 4 (°) 

1 

0 ^ 

£ 


“2. 


h , + Aph s (0^,(0) 


The solution is 




a 

b' 


M 

tt 2. 


c 

d 


1 

E 

sz 


where a = - co 2 3/det & 

b = Ap [0)33 ©24 - ©34 ©23]/ dot £2 
c = © 13 / det £2 

d = Ap [-©33 ©14 + ©34 ©i3]/det £2 
det £2 = ©13 ©24 ■ ©14 ©23 


(3.A34) 


(3.A35) 


and all ©y are evaluated at z = 0. 

We can now calculate the time-dependent motion of relief on the surface and 
crust-mantle boundary by using the kinematic constraints 


dh s 

dt 

dh 

“dt 


= v z (0) * £ a, © j2 (0) 
r = v z (H) = £ a. © i2 (H) 


(A36) 


Using (3.A35), we write (3.A36) as a system of first-order differential equations 
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d_ 

dt 



(3.A37) 


where s-t = aco 12 (0) + c 0)22(0) 

m-i =* bo) 12 (0) + dco 22 (0) - Ap 0)32(0) 
s 2 ss a 0)1 2 (H) + c cd 2 2 (H) 
m 2 = b o) 12 (H) + d co2 2 (H). 


The eigenvalues of the 2 x 2 matrix in (3.A37) are 

1 r 22 i/2i 

X 1 2 = J [ S 1 + m 2 ± (s 1 + m 2 - 2s 1 m 2 + 4s 2 J (3.A38) 

Corresponding eigenvectors are 

x, » (Ci. i) T 1-1.2 


where 



X,-m 


2 


(3.A39) 


The time-dependent solution to coupled topographic decay on the surface and 
crust-mantle boundary may be written 



-tl't 1 -tk 2 

a 1 x 1 e + ag x 2 e 


(3.A40) 



74 


where xj = - 1/Xj . From the initial conditions, the coefficients a-) and a 2 in (3.A40) 
are given by 


a i - h mo - a 2 

(3.A41 ) 

a 2 = (^so ’ Cl hmo)/(C2 ' Cl) 

where h so and h^ are the initial relief at the surface and crust-mantle boundary, 
respectively. The time constants Xj may be expressed in dimensional form by multi- 
plying by Ti 0 /(p c gH), and the amplitude coefficients aj may be redimensionalized by 
multiplying by H. 
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APPENDIX 3.B: 

Support of Crater Relief by Lithospheric Strength 

In Appendix 3.A, the interior of Venus is taken to be a layered viscous fluid 
without finite elastic strength. A simple and commonly adopted alternative 
representation for the mechanical properties of a planetary interior is a uniform thin 
elastic plate overlying an inviscid fluid. With this representation a portion of the 
topographic relief of an impact crater can be supported by bending moments and 
shear stresses within the plate and by buoyancy forces due to displacement of the 
fluid substrate. A lower bound on flexural rigidity or, equivalently, elastic plate 
thickness may be obtained by requiring the flexurally supported relief to be 
consistent with the observed crater depths. The bending stress predicted by this 
elastic plate model, however, exceeds the probable strength in the upper and 
lower portions of the lithosphere. Taking this strength into account, we employ an 
elastic-plastic plate with the same bending moment as the elastic plate to convert 
the lower bound on plate thickness to an upper bound on lithospheric thermal 
gradient. This bound on thermal gradient may then be compared with those 
derived from viscous relaxation models. 

Elastic plate model 

If the topographic relief is small compared with the dominant horizontal 
wavelengths and with the thickness L of the elastic plate, then a crater topographic 
profile h(r) can be considered as a surface load on the plate and the resulting 
vertical deflection w(r) satisfies 


DV4 w (r) + p m gw(r).pgh(r) 


(3.B1) 
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where D is the flexural rigidity, p m is the density of the fluid substrate (taken equal 
to the density of the mantle), p is the density at the surface, g is the gravitational 
acceleration, and w is positive downward. The flexural rigidity is related to the 
plate thickness L by D = E L 3 /[12(1- v 2 )], where E and v are Young’s modulus and 
Poisson’s ratio in the plate, respectively. The final surface topography is given by 
h - w. 

For a given initial crater profile, we seek a minimum flexural rigidity such that the 
crater relief after lithospheric flexure remains consistent with observed crater 
depths. We utilize the same statistical approach to testing the plausibility of elastic 
plate models as for viscous models. Specifically we reject any model that predicts 
a crater depth shallower than the 95% confidence lower bound on mean observed 
crater depth. This criterion is equivalent to a change in relief of no more than 
30-35% for craters with diameters between 30 and 140 km. 

We solve (3.B1), following Solomon and Head [1979], by representing the 
topography as a set of concentric cylinders and utilizing the analytical solution for w 
for a uniform cylindrical load. Because the resulting solution is for the flexure of a 
uniform elastic plate, we must take crustal thickness H to be effectively zero (H « 
L) or infinite (H » L); the corresponding values for p are 3.4 and 2.9 Mg/m 3 , 
respectively. We solve 3.B1) with both E = 10 10 and 10 11 Pa; this range for 
Young's modulus probably spans that of the mean value for the Venus lithosphere, 
with the lower value more likely [Solomon and Head, 1 984], Poisson's ratio for the 
elastic plate is taken to be 0.25. 

Minimum values for elastic plate thickness L for different combinations of crater 
size, lithospheric elastic modulus, and crustal thickness are given in Table 3.B1. 
The areal mean M of the bending moment beneath the crater is also given. An 
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example of the effect of the maximum permitted flexural deformation at 95% 
confidence on a crater topographic profile is shown in Figure 3.B1. 

Elastic-plastic plate model 

Laboratory measurements of the strength of rock suggest that a uniform elastic 
plate is generally a poor model for the mechanical lithosphere. Crustal and mantle 
materials display a finite strength that varies with depth; this strength is controlled 
by frictional sliding on faults at shallow depths and by ductile flow at greater depths 
[Goetze and Evans, 1979; Brace and Kohlstedt, 1980]. In the ductile regime, an 
approximate strength may be estimated by solving equation (3.6) at a specified 
strain rate. The strength depends strongly on temperature (and hence on thermal 
gradient) but only weakly on strain rate. We adopt e = 3 x 10" 17 s -1 (corresponding 
to unit strain in 10 9 yr); a tenfold increase in e results in a change of only a few 
K/km for dT/dz less than 40 K/km. The flow laws for diabase and olivine cited in the 
text are assumed to govern ductile flow in the Venus crust and mantle, respectively. 

In a purely elastic plate, the bending stress varies linearly with depth and is 
antisymmetric about the neutral surface at a depth L/2. For the elastic-plastic plate 
the stress distribution will be modified where the predicted elastic stress exceeds 
the local strength. From the lower bound on elastic plate thickness we obtain an 
upper bound on thermal gradient dT/dz with the elastic-plastic plate model under 
the assumptions that the bending moment of the elastic : plastic plate equals the 
moment M for the purely elastic plate and that the vertical integral of bending stress 
is zero. These conditions lead to two equations that may be solved numerically for 
dT/dz and for the depth z n of the neutral surface: 
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M 



0 


(3.B2) 


J adz = 0 (3.B3) 

o 

where a is the lesser of the strength and the radial horizontal bending stress given 
by the elastic plate model. 

An example of the flexural stress distribution for the elastic-plastic plate 
corresponding to the model in Figure 3.B1 is shown in Figure 3.B2. The upper 
bounds on thermal gradient are given in Table 3.B1. 
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TABLE 3.1 . Large impact craters on Venus with measured depths [Basilevsky et al., 1987] 


Name 

Latitude, 

on 

Longitude, 

°E 

Diameter, 

km 

Depth, 

km 

Elevation a, 
km 

Geological 

Unit b 

Crater 
Class c 

Klenova 

78 

104 

144 

0.78 d 

N/A 

rp 

2 

Cochran 

52 

143 

104 

0.7 

- 0.4 ± 0.4 

sp 

2 

Zhilova 

66.5 

126 

60 

0.62 

0.7 ±0.3 

rp 

2 

LagerlOf 

81 

289.5 

58 

0.62 

N/A 

rp 

3 

Fedosova 

45 

172 

56 

0.75 

- 1 .5 ± 0.3 

sp 

3 

Voynich 

35.5 

56 

56 

0.44 

- 0.5 ± 0.3 

sp 

1 

Yablochkina 

48 

194.5 

55.® 

0.7 

- 0.4 ± 0.2 

rp 

2 

Cotton 

71 

300 

52 

0.7 d 

1.2 ±0.5 

b 

1 

Fedorets 

59.5 

65.5 

52 

0.55 d 

-0.4 ±0.3 

sp 

1 

Duncan 

68 

292 

46 

0.44 

-0.4 ±0.2 

rp 

1 

Tsvetaeva 

65 

147 

45.® 

0.9 

- 0.8 ± 0.3 

sp 

1 

Obukhova 

71 

290 

44 

0.53 d 

0.2 ± 0.3 

b 

2 

Ruslanova 

84 

15 

44 

0.4 

N/A 

rp 

2 

Akhmatova 

61 

308 

42.e 

0.4 

- 0.2 ± 0.4 

rp 

3 

Deledda 

76 

128 

32 

0.5 

N/A 

rp 

2 


a Mean and standard deviation of the elevation measurements obtained by the Pioneer Venus altimeter 
[Pettengill et al., 1980; Ford, 1986] within a circle 2® in diameter centered on the crater, excluding the 
crater interior. The elevation datum is a sphere of radius 6051 .9 km. N/A signifies no data available. 


b "Terrain type" of Barsukov et al. [1986] and Basilevsky et al. [1987]; sp - smooth plains, rp - rolling 
plains, b - ridge and groove belts. 
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c Morphological class [Basilevsky et ai, 1987]: class 1 craters, inferred to be the freshest, have a 
prominent rim and a radar-bright halo; class 2 craters have morphological features similar to class 1 but lack 
a radar-bright halo; for class 3 craters all or most diagnostic morphological features are present but 
degraded. 

d Mean value where two measurements were reported. 


e Geometric mean of two principal diameters. 



TABLE 3.2. Adopted model parameters 


Crater diameter 32 km 

50 km 

140 km 

Wavenumber band 0.060-1 .1 37 km-1 

0.042-1.215 km-l 

0.011-0.325 km-1 

Number of wavenumbers 20 

30 

30 

Average stress difference o 2 MPa 

2 MPa 

4 MPa 

Crustal density 

2.9 Mg/m3 


Mantle density 

3.4 Mg/m3 


Surface temperature 

738 K 


Depth of brittle-ductile transition, crust 

1 km 


Depth of brittle-ductile transition, mantle 

10 km 



Stress evaluated under the assumption that crater depth scales as g- 1 for a reference model 
with crustal thickness H ■ 13 km and mean lithospheric thermal gradient dT/dz - 13 K/km. 
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TABLE 3.B1. 

Limits on plate thickness and thermal gradient beneath Venus impact craters 
from elastic and elastic-plastic models for lithospheric flexure 



Elastic Plate 


Elastic-Plastic Plate 

D = 140km 

H 

E, 

L, 

M, 

dT/dz, 


Pa 

km 

10" N 

K/km 

0 

ioio 

20 

4.2 

19 


10" 

9 

9.0 

29 

H »L 

10l° 

22 

4.8 

5 


10" 

10 

10.3 

6 

D = 50km 

H 

E, 

L, 

M, 

dT/dz, 


Pa 

km 

10" N 

K/km 

0 

10l° 

7.7 

2.0 

49 


10" 

3.5 

1.9 

68 


H»L 

ioi° 

7.9 

2.1 


10" 

3.7 

2.1 


18 
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D » 32 km 


H 

E, 

U 

M, 

dT/dz, 


Pa 

km 

1011 N 

K/km 

0 

10 10 

3.3 

1.1 

96 


ion 

1.5 

2.3 

110 

H»L 

ioio 

3.5 

1.2 

31 


ion 

1.6 

2.6 

30 
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FIGURE CAPTIONS 

Figure 3.1 Venera 16 radar image of Klenova, an impact structure 144 km in 
diameter [Barsukov et al., 1986]. Klenova has a complex morphology, 
including a scalloped rim crest, a flat inner floor, and a peak ring. The 
width of the image is 520 km; white lines denote meridians spaced at 
5° intervals; the 80° parallel is also shown. 

Figure 3.2 Depth versus diameter for impact craters on Venus, from Basilevsky et 
al. [1987]. Solid circles denote craters included in this study (Table 
3.1); crosses denote other craters not included here, for reasons 
given in the text. The cross at upper right represents Cleopatra 
Patera. Also shown are the mean and 95% confidence limits for the 
mean depths of observed (solid lines) and predicted fresh (dashed 
lines) impact craters. Depths of fresh craters are estimated from the 
depths of lunar Copernican and Eratosthenian craters under the 
assumption that relief scales inversely with surface gravitational 
acceleration (see text). 

Figure 3.3 Locations of the 15 impact craters used in this study (filled circles). 

Topographic provinces, following Masursky et al. [1980], are outlined; 
highlands are cross-hatched and lowlands are shaded. Most of the 
craters depicted lie at intermediate elevations. 

Figure 3.4 Distribution of impact craters with depths measured by Venera 1 5 and 
16 (from Table 3.1). (a) Distribution by regional elevation as indicated 
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by Pioneer Venus altimetric data (where available), (b) Distribution 
by geological unit [Basilevsky et ai, 1987], where sp denotes smooth 
plains, rp rolling plains, and b ridge and groove belts, (c) Distribution 
by crater morphological class [Basilevsky et ai, 1987], where class 1 
craters are most pristine and class 3 are significantly degraded. The 
craters used in this study are fairly typical in terms of elevation, 
geological terrain, and state of preservation. 

Figure 3.5 The statistical approach to viscous relaxation adopted in this paper. 

Assuming that at a fixed crater diameter D the crater depths follow a 
log-normal distribution, and that all craters at diameter D have an age 
given by half the mean crater retention age, then a particular viscous 
relaxation scenario will reduce the mean depth of fresh craters from d-, 
to 5. If 5 is much smaller than d v , the observed mean depth of craters 
of diameter D, then that viscous relaxation model is not plausible. A t 
test may be used to quantify this relationship by calculating the 
significance level a at which the null hypothesis d v * 8 can be 
rejected in favor of d v > 8; larger values of a are associated with more 
plausible viscous relaxation models. 

Figure 3.6 Two models for the viscous relaxation of a Venus impact crater 140 
km in diameter, illustrating the combined effects of thermal gradient 
(dT/dz) and crustal thickness (H) on the extent of relaxation versus 
crater age t. The initial topography is scaled from that of the lunar 
crater Langrenus and is assumed to be initially uncompensated (c = 
0). (a) A model with dT/dz = 10 K/km and H = 10 km. There is little 
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Figure 3.7 


i 


) 

Figure 3.8 


I 


change in relief after 50 m.y., but some decrease has occurred by 500 
m.y. The latter change could be marginally resolved at 95% 
confidence (see text), (b) A model with dT/dz = 1 5 K/km and H = 1 5 
km. Substantial topographic relaxation occurs in 50 m.y., and the 
crater relief is effectively removed by 500 m.y. 

Two models for the viscous relaxation of a Venus impact crater 32 km 
in diameter. Except for crater diameter and initial topography (scaled 
from that of the lunar crater Timocharis), parameters are identical to 
those of the models depicted in Figure 3.6. (a) A model with dT/dz * 
10 K/km and H = 10 km. The change in relief is small, even after 500 
m.y. (b) A model with dT/dz = 15 K/km and H = 15 km. The extent of 
viscous relaxation is greater than in (a) but the change in relief is not 
proportionally nearly as great as for the crater with D = 140 km, in part 
because the relaxation of larger features is sensitive to the generally 
lower viscosity at greater depth. 

Relaxation time constants x versus wavenumber k for the models in 
Figures 3.6 and 3.7. The solid portions of each curve delimit the 
wavenumber range corresponding to the largest values of the crater 
topographic spectra (2rc/k = 0.5-3 D). At dT/dz = 1 5 K/km and H = 15 
km, the relevant time constants for D = 32 km exceed those for D = 
140 km. At dT/dz = 10 K/km and H = 10 km, however, the relevant 
time constants for both diameters are comparable because the 
sensitivity of the smaller crater to the viscosity of the weak lower crust 
is largely offset by an increase in effective viscosity due to lower 
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Figure 3.13 Summary of 95% confidence levels for rejecting combinations of 
dT/dz and H from Figures 3.9-12 (solid lines) for craters 50 m.y. (a) or 
500 m.y. (b) in age. A similar curve is shown for a crater 50 km in 
diameter; the initial topography is scaled from that of the lunar crater 
Eratosthenes. Dashed extensions of several of these curves show 
similar calculations for initially compensated craters where these 
solutions differ from those for initially uncompensated craters; the 
characteristic stress difference o and thus the effective viscosity at 
each temperature are taken to be identical to those for the initially 
uncompensated cases. The dashed curve at upper right gives the 
95% confidence level for rejecting combinations of dT/dz and H for an 
initially compensated crater 140 km in diameter and a characteristic 
stress difference of 0.5 MPa. Such a low value of o is unlikely but 
gives conservative upper limits to dT/dz and H for initially 
compensated large craters. 

Figure 3.A1 . Sketch of the viscous relaxation problem as applied to Venus in this 
paper. A crust with density p c and thickness H overlies a mantle with 
density p c + Ap. Given an average lithospheric thermal gradient 
dT/dz and specified flow laws for ductile deformation of crustal and 
mantle rocks, the effective viscosity as a function of depth may be 
constructed. 

Figure 3.B1. Flexural deformation of a 140-km-diameter impact crater on an elastic 
plate overlying an inviscid fluid. The portion of the original crater 
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differential stress. 


Figure 3.9 Effect of viscous relaxation for an initially uncompensated crater 140 
km in diameter after 50 m.y., shown as a function of dT/dz and H. (a) 
Fractional relief remaining, (b) Confidence levels for rejecting 
combinations of dT/dz and H. At a specified confidence level, models 
to the right of each curve may be rejected. 

Figure 3.10 Effect of viscous relaxation for an initially uncompensated crater 140 
km in diameter after 500 m.y., shown as a function of dT/dz and H. (a) 
Fractional relief remaining. The small shift in the 0.1 contour is due to 
a change in the radial position of the point of greatest depth on the 
crater floor, (b) Confidence levels for rejecting combinations of dT/dz 
and H. 

Figure 3.11 Effect of viscous relaxation for an initially uncompensated crater 32 
km in diameter after 50 m.y., shown as a function of dT/dz and H. (a) 
Fractional relief remaining, (b) Confidence levels for rejecting 
combinations of dT/dz and H. A greater proportional relief remains 
than under the same conditions with D = 140 km (Figure 3.9). 


Figure 3.12 Effect of viscous relaxation for an initially uncompensated crater 32 
km in diameter after 500 m.y., shown as a function of dT/dz and H. (a) 
Fractional relief remaining, (b) Confidence levels for rejecting 
combinations of dT/dz and H. Note that the extent of relaxation is 
insensitive to crustal thickness for H greater than about 20 km. 
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relief (solid line) that can be elastically supported increases with the 
plate thickness. The dashed profile corresponding to an elastic plate 
thickness of 20 km for E = 10 10 Pa matches the 95% confidence lower 
bound on the mean observed depth of craters of this size. 

Figure 3.B2. The distribution with depth of the average radial horizontal stress 
beneath the crater in the example shown in Figure 3.B1 for an 
elastic-plastic plate model. The crustal thickness is taken to be zero, 
so the bending stresses are limited by the strength envelope for 
olivine (long dashed line). An upper bound on the thermal gradient 
for this model is 19 K/km. 
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INTRODUCTION 

Aphrodite Terra plays a key role in the geology of Venus. It is the largest 
highland region on the planet, stretching nearly 15,000 km in length [Pettengill et 
al., 1980]. Aphrodite Terra together with Beta Regio and its extensions consititute 
the equatorial highlands, which stand out on the basis of their rough, elevated 
topography, large free air gravity anomalies, and volcanic and rift structures 
inferred from radar observations [Masursky et at., 1980; Phillips et al., 1981]. 
Recently-mapped features in Aphrodite include an organized pattern of lineaments 
[Crumpler et al., 1987] and elements of topographic bilateral symmetry [Crumpler 
and Head, 1988a] in a direction parallel to these lineaments. 

Upwelling of mantle convective flow is a likely explanation for most of these 
observations, but the planform of this convection and its interaction with the surface 
is unresolved. Many workers [Phillips and Malin, 1983; Morgan and Phillips ; 1983; 
Banerdt, 1986; Kiefer et a!., 1986] have proposed mantle-plume and hotspot 
models, from inferences that topography is maintained by mantle convection and 
that regional tectonics are largely vertical with little lateral strain. Head and 
Crumpler [1987], however, have suggested that Aphrodite is a divergent plate 
boundary, on the basis of geological evidence for horizontal displacements over 
several thousand kilometers. 

In this paper, we first review the characteristics of Aphrodite Terra and their 
relation to these tectonic models. We suggest an additional model in which linear 
upwelling and divergence of mantle flow is accommodated at the surface by 
distributed, nonrigid deformation. We undertake several quantitative tests of 
evidence proposed to support the rigid-plate model, namely, that the lineaments, 
inferred to be fracture zone traces, are organized in a geometrical pattern that may 
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be interpreted in terms of plate kinematics, that the topography is consistent with a 
diverging, cooling thermal boundary layer, and that a regional topographic bilateral 
symmetry exists. Since terrestrial analogy is an important part of the 
plate-divergence hypothesis, we perform the same tests on the Mid-Atlantic Ridge, 
restricting our information to a level comparable to that for Venus. We examine the 
implications of these tests for the other models as well, and we suggest 
observations from the planned Magellan mission that will lead to the resolution of 
these competing hypotheses. 

CHARACTERISTICS OF APHRODITE TERRA 

Regional elevations exceeding 2 km above the mean planetary radius are 
restricted to the three principal highlands of Venus: Ishtar Terra, Beta Regio, and 
Aphrodite Terra [Masursky et ai, 1980]. Aphrodite and Beta may be distinguished 
from Ishtar on the basis of their smoother overall topography, larger free-air gravity 
anomalies, and a predominance of tectonic features of inferred extensional origin 
(see below). 

Aphrodite Terra consists of three mountainous plateaus (Ovda, Thetis, and Atla 
Regiones) and intervening equatorial terrain between longitudes 50° and 210° 
(Figure 4.1). Thetis and Atla are connected by a broad saddle which is dissected 
by numerous troughs and ridges; we refer to this region as east-central Aphrodite. 
The principal tectonic elements of Aphrodite are linear to arcuate troughs, or 
chasmata, often flanked by asymmetric ridges, which broadly resemble terrestrial 
rift valleys. Groups of aligned paired troughs and ridges become clearly resolvable 
in Pioneer Venus (PV) altimetry over length scales of several hundred kilometers, 
and arcuate troughs and ridges range up to several thousand kilometers in length 
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and several kilometers in depth [McGill et al., 1983]. The largest of these features 
occur on the southern flanks of Ovda and Thetis Regiones and in east-central 
Aphrodite. The single largest trough, Artemis Chasma, is perhaps the most 
anomalous, in that it is strongly curved, extending over 210° of arc. Chasmata are 
also observed in the Beta-Phoebe region; detailed analysis of this area from 
Pioneer Venus [McGill et al., 1981], Arecibo [Campbell et al., 1984] and limited 
Venera 15/16 data [Stofan et al., 1987] suggests that these features are graben 
and therefore the result of extensional stresses. 

In addition to the ridge-and-trough systems, two other classes of potential 
tectonic features that may be related to the evolution of the equatorial highlands 
have been described by Masursky et al. [1980]. The first includes lineaments and 
scarps associated with highlands that may indicate terraced block faulting, such as 
a series of radar-bright and dark bands on the northwest margin of Ovda. The 
second consists of linear "disturbed zones" of modest topographic relief (several 
hundred meters) that are rough at radar wavelengths; such zones link Atla Regio to 
Beta Regio and Themis Regio. The association of these disturbed zones with 
troughs lead Schaber [1982] to infer an extensional origin for them as well. 
Schaber further proposed that elevated regions (Beta, Phoebe, Atla) are located at 
sites of weakened lithosphere at the intersections of major zones of lithospheric 
disruption. 

Long-wavelength gravity and topography are correlated on Venus [Phillips et 
al., 1979]. This correlation is strongest for Atla and Beta Regiones, and weaker for 
Ovda and Thetis Regiones [Sjogren et al., 1983]. Indeed, there appears to be a 
progressive variation in the gravity signature from Ovda to Atla, manifested by an 
increasing apparent depth of isostatic compensation across this region [Herrick 
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and Bills, 1988]. 

Phillips et al. [1981] analyzed the trade-off between flexural and crustal (Airy) 
compensation in western Aphrodite and found that static support of relief requires 
either an elastic lithosphere several hundred kilometers thick or a crustal thickness 
compensation depth of about 100 km. If heat production in the Venus interior is 
comparable to that of Earth, high temperatures at these depths lead to solid-state 
creep and relaxation of stress and topography for both cases. Phillips et al. [1981] 
concluded either that the crust and upper mantle of Venus are depleted in 
heat-producing elements or that the topography is young or maintained by interior 
convective flow. 

Others have sought to explain the observed topography by density variations 
within the mantle. Morgan and Phillips [1983] found that 99% of the topography of 
Venus could be explained by a hotspot model in which relief is primarily generated 
by thermal uplift due to lithospheric thinning, supplemented by volcanic 
construction. Banerdt [1986] presented a model that matches topography by 
variations in upper mantle density and crustal thickness and which explicitly 
satisfies the observed gravity field. He concluded that uplift by low-density mantle 
is generally required to produce tensional stresses normal to the strike of troughs in 
Aphrodite and the disturbed zones joining Aphrodite to Beta and Themis Regiones. 
Banerdt [1988] also argued that these density differences are too deep to be 
statically maintained and therefore must be caused by temperature variations 
associated with mantle convection. Within Aphrodite itself, a thicker crust and 
lesser mantle density anomaly (presumably the result of a cooler upper mantle) are 
required under this model for Ovda Regio than elsewhere. 

The dynamic support of topography (i.e., including a component of support by 
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normal stresses induced by upward mantle flow) by mantle plumes on Venus has 
been investigated by Kiefer etal. [1986] and Kiefer and Hager [ 1987]. They found 
that models without an upper-mantle low-viscosity zone allowed stronger coupling 
of convective stresses to the lithosphere, producing geoid anomalies and topo- 
graphy comparable to those observed in the equatorial highlands. Kiefer and 
Hager [1987] further suggested that the comparatively smaller geoid anomaly at 
Ovda Regio is due to a hotter upper mantle and a better-developed low-viscosity 
zone, leading to a smaller admittance and apparent depth of compensation. 
Robinson et al. [1986] described a similar effect for terrestrial mid-plate swells, 
wherein an upper-mantle low-viscosity zone could explain apparent depths of 
compensation within the lithosphere. These analyses suggest both that the mantle 
actively contributes to topographic support and that lateral variations in the style of 
support exist across Aphrodite. 

Recent detailed geological mapping of Aphrodite Terra has suggested two new 
classes of tectonic features. From Pioneer Venus orbital radar altimetric and 
imaging maps and from Arecibo radar altimetric profiles, Crumpler et al. [1987] 
mapped a series of parallel linear structural discontinuities in Ovda and Thetis 
Regiones, termed "cross-strike discontinuities", or CSDs, which are up to several 
thousand kilometers in length and show an approximately constant NNW strike, at 
high angles to the general trend of Aphrodite. Mapping of additional CSDs, 
extending through east-central Aphrodite, was reported by Crumpler and Head 
[1988b]. A strong bilateral symmetry across Aphrodite Terra parallel to the strike of 
the CSDs was noted by Crumpler and Head [1 988a], not only in the broad shape of 
the highland but also in individual topographic elements to the north and south. 
Assembling the centers of symmetry over many profiles, Crumpler and Head 
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argued that these centers were organized along a series of linear, offset segments. 
The system of CSDs and linear centers of topographic bilateral symmetry as 
suggested by these authors is shown in Figure 4.2. 

TECTONIC MODELS FOR APHRODITE TERRA 

We describe three models for the tectonics of Aphrodite Terra and its interaction 
with mantle convection. Vertical tectonics will dominate if little horizontal surface 
displacement results from an upwelling mantle flow. Rigid-plate divergence is the 
appropriate description if shear strain from large horizontal displacements is 
accommodated only in narrow deformation zones. Alternatively, such strain could 
be broadly accommodated if the lithosphere is thin and easily deformed. We refer 
to these models as vertical tectonics, plate divergence, and distributed deformation, 
respectively (Figure 4.3). Specific predictions of these models are given below. 

Vertical Tectonics 

If divergence of upward mantle flow is not strongly coupled to the surface, 
predominantly vertical tectonics will result. Individual mantle plumes beneath 
Ovda, Thetis, and Atla Regiones can account for their approximately circular 
topographic plan [Kiefer et al., 1986], and no link to global tectonics is required. 
Axisymmetric plumes do not adequately account for the east-west orientation of 
chasmata, however, unless individual plumes are contemporaneous, so that a 
regionally coherent north-south orientation of extensional stress may develop. In 
this situation, there may be little difference between a chain of simultaneous, 
linearly aligned plumes and a linear zone of upwelling appropriate to the 
alternative models. 
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The consistent orientation of the CSDs (nearly normal to the trend of Aphrodite) 
also cannot be explained by widely separated plumes. Even for a linear upwelling, 
the great lateral extent of the CSDs implies that horizontal strains exist to distances 
of several thousand kilometers from the topographic axis of Aphrodite Terra. 
Although it is possible that topographically-induced stresses may be transmitted to 
such distances [e.g, Bindschadler and Parmentier, 1987], such stresses should 
produce compressional features oriented approximately parallel to the highland 
trend and not nearly perpendicular features, as observed for the CSDs. 

If upwelling occurs along linear zones, then the divergence and cooling of new 
lithosphere can lead to subsidence of topography. Because the surface is 
decoupled from the interior in the vertical-tectonics model, however, the thermal 
structure will differ from that in cases where the surface is part of a plate, e.g., 
terrestrial oceanic lithosphere. This in turn will lead to large departures from the 
form of topography expected from a thermal boundary layer (see below). 

Plate Divergence 

The CSDs are the principal geological indicator of large-scale horizontal motion 
in Aphrodite. Crumpler et al. [1987] found that the pattern of CSDs was consistent 
with either regional strike-slip faults or oceanic fracture zones, but they favored the 
latter interpretation on the basis of the parallelism and great length of the CSDs. 
Crumpler and Head [1988a] mapped individual elements of topographic bilateral 
symmetry along profiles parallel to the CSDs, and they suggested that these 
smaller topographic features were created at the divergent plate boundary and 
subsequently rifted and separated. 

As additional support of the plate-divergence hypothesis, Head and Crumpler 
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[1987] argued that the long-wavelength relief of Aphrodite Terra is consistent with 
divergence, cooling and subsidence of newly created lithosphere. For a constant 
rate of divergence, topographic subsidence away from the rise crest should be 
proportional to the square root of distance x[e.g., Parsons and Sclater, 1977]. 
Evidence for such thermal boundary layer topography on Venus was examined by 
Kaula and Phillips [1981]. On a broader scale, they found that the equatorial 
highlands satisfied two out of four important criteria for plate tectonics; namely, that 
the rises were predominantly concave upward and rose above a broad reference 
plain of little topographic variation. However, Kaula and Phillips concluded that 
these highlands do not have a well-defined maximum elevation, nor do they form a 
globally interconnected system, in contrast to Earth's mid-ocean ridges. The first of 
these results may not be a valid argument against the plate-divergence hypothesis 
if there are variable contributions to topography from dynamic support or crustal 
thickness variations. The second argument, however, remains a serious objection 
to the model, as other plate boundaries integrating Aphrodite into a global system 
have not been identified. On the basis of thermal boundary-layer scaling 
arguments, Kaula and Phillips [1981] predicted that topographic relief at divergent 
margins on Venus should be about 40% that of Earth's mid-ocean ridge elevations. 
The equatorial highlands have relief in excess of 4 km with respect to the mean 
elevation, in disagreement with this prediction by a factor of two or more. In 
addition to the possible variations in relief just mentioned, this result depends on 
the proportionate scaling between Earth and Venus of heat sources that limit the 
thickness of the boundary layer and amount of subsidence. Greater relative 
cooling and thickening of the thermal boundary layer on Venus, and greater relief 
than that predicted by Kaula and Phillips, cannot be ruled out. We return to these 
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issues in the discussion of thermal boundary layer topography below. 

The chasmata do not play an important role in the plate-tectonic model outlined 
by Crumplerand Head [1988a], in spite of the great size and geographic extent of 
these features. The linear centers of crustal divergence defined by these workers 
generally do not coincide with the principal chasmata. However, this is not 
necessarily evidence against the model, as the orientation of the chasmata is 
generally perpendicular to the inferred direction of plate divergence, and horizontal 
extensional stresses may be expected for this high-standing region. 

Distributed Deformation 

Phillips [1986] analyzed the coupling of convective tractions to a thin 
lithosphere and concluded that such stresses could reach 1 kbar and lead to 
tectonic deformation. Under the distributed-deformation model, then, CSDs may 
represent flowlines indicating the planform of underlying mantle convection. As in 
the plate model, the east-west orientations of the chasmata are consistent with a 
linear axis of upwelling and divergence, but topographic stresses may also 
contribute to their formation. A thermal boundary-layer component of topography 
should also be present, although the apparent rate of spreading would apply 
largely to the underlying mantle and not necessarily to the spreading rate, if any, of 
the surface. Moreover, near-surface thermal structure differing from a simple 
cooling boundary layer could also lead to departures from the Vx subsidence of 
topography. Because surface deformation is not required to be spatially coherent, 
the global interconnection of tectonic features need not be as strongly developed 
as under the plate-tectonic model. 
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TESTS OF TECTONIC MODELS 

If the CSDs are analogs of terrestrial fracture zones, then they provide strong 
geometrical constraints on plate motion. The portions of CSDs offsetting adjacent 
linear centers of bilateral symmetry represent transform faults under the 
plate-tectonic hypothesis and must follow small circles about an appropriate 
rotation pole, and the distal traces of CSDs, representing inactive fracture zones, 
can yield information on past motions. If more than two plates exist, then one or 
more triple junctions must also be present, and additional plate boundaries must 
be identified. Because there is no evidence to support a ridge or trench 
intersection with Aphrodite proper, the most likely such boundary would be a fault 
zone, presumably one of the mapped CSDs. Transform motion would probably 
dominate, although some extension or compression most likely would also occur 
across the boundary. Under nonrigid divergence, CSDs need not satisfy the same 
criteria as do fracture zones. Our first test then, is to examine the CSDs from the 
perspective of plate kinematics. 

Provided that thermal contributions to isostasy dominate at wavelengths of less 
than a few thousand kilometers, the plate-divergence model should show relief that 
falls off approximately with the square root of age. Furthermore, rigid-plate 
rotations constrain these velocities to vary regularly with position along the 
plate-boundary distance. 

The thermal boundary-layer tests of Kaula and Phillips [1981] centered on the 
differences in the distribution of apparent rise heights on Earth and Venus; they 
calculated apparent spreading rates as a function of the deviation of individual rise 
crests from the modal rise height under varying definitions. The variation of 
apparent spreading rate with distance from the rise crest and under conditions of 
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restricted uncertainty in the topography were also tabulated. As mentioned above, 
however, and discussed in detail below, rise height may not be a particularly 
diagnostic guide to the thermal boundary layer contribution to topography. Given 
the specific plate-divergence scenario proffered by Head and Crumpler [1987], it is 
appropriate to extend the analysis of Kaula and Phillips. Our second test is to 
perform least-squares fits of topography versus the square root of distance along 
the spreading direction indicated by Crumpler and Head[ 1988a]. We test the 
goodness of fit by geographic location, both as a function of distance away from the 
rise axis and position along the axis. These tests follow from suggestions by 
Crumpler and Head that regions away from the Ovda and Thetis plateaus, such as 
the rolling plains to the north and south or the saddle-shaped region of east-central 
Aphrodite, may show a stronger thermal boundary layer signature. 

Lastly, we test the bilateral symmetry of topography at short wavelengths by a 
cross-correlation technique, as such relief could support a plate-divergence model 
[Crumpler and Head, 1988a]. On Earth, magnetic lineations have been the primary 
signature of symmetric spreading [e.g., Vine, 1966], although a preliminary catalog 
of potentially rifted and separated topography has been compiled by Crumpler and 
Head [1988c]. It is important, therefore, to examine this question for the terrestrial 
seafloor as well. 

PLATE KINEMATIC TESTS 

By Euler's theorem, the relative motion of two rigid plates on the surface of a 
sphere may be described by a rotation about an axis passing through the center of 
the sphere [Bullard et al, 1965]. The intersection of the axis of instantaneous 
relative angular velocity with the planetary surface is known as the Euler pole. 
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Because transform faults follow small circles about this pole, a single best-fitting 
Euler pole may be calculated by minimizing % 2 [LePichon, 1968]: 



( 4 . 1 ) 


where Ci and q 2 are the azimuth and its uncertainty for the i th transform fault, and 
Coi is the azimuth of a small circle about an Euler pole at latitude Xq and longitude 
(t> 0 . In practice, a search over all (Xq, <j>o) is used to find the minimum value of x 2 , 
and formal errors corresponding to different values of % 2 may be assigned 
following standard statistical methods [e.g., Press etal., 1986]. 

If the relative position of the Euler pole and the plate boundary does not change 
with time, then fracture zones should follow the same small circles as do transform 
faults. On Earth, however, it is well documented that Euler poles do migrate 
relative to the plate boundary over time scales of millions of years. We therefore 
compute an instantaneous Euler pole using only the segments of CSDs located 
between offset linear centers of bilateral symmetry (Table 4.1). The calculated 
Euler pole lies at 20 ±3°N, 204 ± 17°E. Because of the fewer number of CSDs in 
east-central Aphrodite and their larger location errors, the Euler pole position is 
fixed largely by the CSDs in Ovda and Thetis Regiones. Figure 4.4 shows the 
location of this Euler pole in relation to the Aphrodite system of CSDs and 
symmetry centers, along with level curves for one and two standard deviations. 
Small circles about the Euler pole are concave to the northeast in eastern 
Aphrodite and to the southwest in western Aphrodite. This is a consequence of a 
pole location approximately 90° from the center of the CSD system, which in turn is 
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due to the fact that the CSD traces are all nearly parallel in Mercator projection. 

The geometry of the CSDs in relation to the best-fitting Euler pole is further 
elucidated by examining CSD geometry in a spherical coordinate system about the 
Euler pole (Figure 4.5). Small circles about the Euler pole now appear as lines of 
latitude in this projection, and lines of longitude are great circles passing through 
the Euler pole. The linear centers of topographic bilateral symmetry fall 
approximately along lines of longitude for Euler latitudes 40°S to 35°N. Within this 
latitude range, CSD segments within about 30° of the axes of topographic 
symmetry tend to follow parallels. The CSDs display increasing curvature, 
however, at larger Euler latitudes. The opposing curvature of CSDs on the same 
side of the axis of topographic symmetry is evidence against a single, fixed Euler 
pole, as movement of rigid lithosphere cannot be geometrically accommodated at 
increasing distances from the rise axis. A single Euler pole fixed with respect to the 
plate boundary over the age of the CSDs may therefore be rejected. We examined 
the effect of uncertainties in CSD strikes on this result by posing a t test for the 
deviations of CSD azimuths from lines of latitude in the Euler-pole projection. For 
CSD segments greater than 20° distant from both the Euler equator and from the 
nearest linear center of topographic bilateral symmetry, the inward curvature may 
be resolved at greater than 95% confidence. This conclusion may not be valid, 
however, if there is some systematic, rather than random, measurement error. 

It is not unreasonable that a single, fixed Euler pole should not match the 
observed CSD traces, simply given their size and probable age. The Aphrodite 
system of linear centers of topographic bilateral symmetry is over 12,000 km in 
length. By comparison, the longest plate boundaries on Earth are less than 10,000 
km long. In particular, x 2 = 0.07 for the best-fitting Euler pole above; a similar test 
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on the Atlantic Ocean basin using 13 transform faults between latitudes 52°N and 
30°S with offsets exceeding 1° yielded a pole at 55 ± 12° N, 328 ± 2° E, and x 2 * 
0.08, assuming measurement errors comparable to those found for the CSDs. In 
reality, the Atlantic system includes four major plates whose speed and direction of 
movement have varied over time. 

The next simplest model is a single rotation pole, but with time-dependent 
motion with respect to the plate boundary. On Earth, for example, migration of the 
pole with respect to the Mid-Atlantic Ridge (MAR) in the late Mesozoic and early 
Cenozoic has imparted a broad sigmoidal shape to fracture zones on either side of 
the MAR [e.g., Klitgord and Schouten, 1986]. If changes in the location of the pole 
for a two-plate system in Aphrodite have occurred, then all fracture zones on one 
side of the ridge system should show a shift in azimuth in the same sense, and the 
shifts on the opposite side should be in the opposite sense. The opposing 
curvature of CSDs on the same side of the axis of topographic symmetry in Figure 
4.5 argues against a single Euler pole, regardless of whether it migrates with 
respect to the plate boundary over time or not. 

The lack of Pioneer Venus radar images further to the east (in Euler 
coordinates) precludes consideration of more complicated multiple-pole, 
time-dependent models. Although such models cannot be tested quantitatively on 
the basis of the data considered here, any scenario would require a separate pole 
for each few CSDs, or else poles in a system of several plates would have to be 
adjusted continuously with time in order to maintain the fairly constant strike of the 
CSDs as portrayed in Mercator projection. The former scenario may probably be 
ruled out because it implies that each CSD is a plate boundary and would lead to 
plates that are very long in the spreading direction compared with the length of the 
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divergent boundary; stresses on such a "ribbon" plate would most likely break it 
into smaller, more equidimensional pieces, leading to more distributed 
deformation. Therefore, if the plate-divergence model holds, the Aphrodite system 
must be limited to only a few plates and Euler poles must move with respect to 
plate boundaries over time to maintain the nearly-linear traces of the CSDs. We 
predict additional features of such a geometry in a later section. 

SELECTION AND PROCESSING OF TOPOGRAPHIC DATA 

Some preliminary processing of topographic data is necessary prior to 
conducting tests for thermal boundary layer topography and bilateral symmetry. 
Each test of topographic characteristics of Aphrodite Terra is matched by an 
equivalent test of the MAR; the latter data are restricted to Pioneer Venus resolution 
by the following approach. There are nearly 20,000 individual PV radar altimetric 
measurements over the greater Aphrodite area ( approximately 130x10 6 km 2 ) 
depicted in Figure 4.1 [Pettengill et ai, 1980; Ford, 1986], corresponding to a mean 
spacing of ~ 80 km. The geometric mean of the along- and cross-track dimensions 
of PV altimeter footprints at typical altitudes over Aphrodite is ~ 25 km. The 
terrestrial DBDB5 data set [U.S. Naval Oceanographic Office, 1985] consists of 
topographic and bathymetric data interpolated every 5’ of arc, which may be 
considered equivalent to an 8 km footprint and spacing at 30° latitude. It is 
appropriate, therefore, to decimate these data by sampling one 1 5’xl 5* square out 
of each square degree to simulate the coverage and resolution of Pioneer Venus. 

We examined the elevations along 38 profiles across Aphrodite and 35 profiles 
across the MAR. For Aphrodite, the traces of these profiles are taken perpendicular 
to the linear centers of topographic bilateral symmetry and follow great-circle paths 
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approximately parallel to the CSDs (Figure 4.6). Profiles centers are spaced 
approximately 200 km apart (Table 4.2) and are taken no closer than about 150 km 
to CSDs. All data within 100 km of the profile trace are used; therefore 
near-contiguous coverage of Aphrodite Terra away from CSDs is afforded. For the 
Atlantic, profile centers were also spaced at approximately 200 km intervals along 
the mid-ocean ridge (Table 4.3) and were taken no closer than about 150 km from 
large-offset (>1°) transform faults. This offset length was chosen on the grounds 
that such transforms might be recognized at a resolution for Earth comparable to 
that available on Venus. The terrestrial profiles approximately follow fracture zone 
traces [Klitgord and Schouten, 1986; Cande at al., 1988] and appear curved 
because of changes in plate motion (Figure 4.7). A few short profiles in the 
equatorial and northern Atlantic follow small circles about the instantaneous 
two-plate rotation pole [Minster and Jordan , 1978]. The mean length of the MAR 
profiles is 1 500 km in each direction from the ridge. 

The Venus profiles are taken to a distance of 3600 km in each direction, 
comparable to those shown by Crumpler and Head[ 1988a] and extending slightly 
beyond the mapped length of the CSDs. Because continents play no role in the 
plate-divergence model for Aphrodite Terra, terrestrial profiles are truncated where 
they intersect continental rises or oceanic plateaus with continental basement. No 
such constraints are imposed on any profiles crossing seamounts or mid-plate 
swells, for which plausible analogs may exist on Venus, regardless of any 
plate-tectonic environment. 

For the analysis of boundary-layer topographic characteristics, all data within 
the 200-km-wide profile swath are taken to correspond to the distance of the 
nearest point on the profile trace. For the cross-correlation approaches used to test 
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bilateral symmetry, data gridded at regular intervals are required. We used a 
100-km grid spacing, comparable to the typical measurement spacing of the PV 
altimeter. For each grid point, elevation measurements within 100 km were 
averaged after weighting by the inverse squared distance. The inverse square 
form is appropriate for the case of sample (footprint) sizes smaller than the 
observation spacing because it strongly favors nearby observations but will include 
distant observations if no nearby observations are available. The topographic 
profiles so constructed for Venus are shown in Figure 4.8, and those for Earth are 
given in Figure 4.9. 

TESTS FOR THERMAL BOUNDARY-LAYER TOPOGRAPHY 

A divergent thermal boundary layer may contribute to relief under any of the 
models discussed here, but its form is best known for the plate-divergence model. 
On Earth, the variation of seafloor depth with age has been attributed to the cooling 
of spreading oceanic lithosphere. Two alternative models exist for the form of this 
cooling and differ in their predictions for the depth of old seafloor. Turcotte and 
Oxburgh [1967] described the thickening of the lithosphere as the thermal 
boundary layer of a mantle convective system. Since this boundary layer is 
formulated in terms of a conductively cooling halfspace, it grows in proportion to the 
square root of distance x until a zone of downwelling is reached. Subsidence is 
also proportional to Vx. In the alternative plate model [McKenzie, 1967; Parsons 
and Sclater, 1977] the thickness of mature lithosphere is limited by some constant 
upper mantle temperature; the plate grows asymptotically to a maximum thickness 
determined by this temperature. Since early cooling resembles that of a halfspace, 
subsidence is initially proportional to Vx, but older seafloor subsides more slowly, 
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decaying exponentially to a constant depth [Parsons and Sclater, 1977]. This 
depth is about 6400 m, and terrestrial seafloor elevations e given herein are 
measured relative to this datum. Where seafloor age is known as a function of 
distance, the variation of seafloor depth may be predicted to well within the 
300-400 m accuracy of regional bathymetric measurements, after correcting for any 
sediments and their isostatic adjustment [Parsons and Sclater, 1 977]. 

We use elevation-distance measurements to infer average apparent spreading 
rates, employing the inner asymptotic solution for the Vx subsidence of diverging, 
cooling lithosphere given by Parsons and Sclater [ 1977]. Where this approx- 
imation is valid, it may apply to either the thermal boundary layer or plate models, 
and so the limits to this form when applied to Venus, if any, must be evaluated. 
There are two potential reasons for deviation from this relation: a limiting plate 
thickness and other sources of topography. 

On Earth, the exact cause of a limiting plate thickness is not completely 
understood. Parsons and McKenzie [1978] interpreted the departure from a Vx 
relation at approximately 70 m.y. in terms of convective instability at the base of the 
lithosphere, leading to small-scale convection and to the supply of additional heat. 
However, Buck [1985] has concluded from temperature-dependent viscosity 
models that small-scale convection may develop much earlier than inferred by 
Parsons and McKenzie from constant-viscosity models. An implication of this result 
is that the flattening of the depth-age curve must have another cause. Other 
mechanisms have been proposed, such as internal heating [Forsyth, 1 975; Jarvis 
and Peltier, 1982] or shear heating [Schubert et ai, 1976]; indeed, Heestand and 
Crough [1981] contend that the depth-age flattening is an artifact of including 
hotspot traces in averages of seafloor depth at large age. 
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Nevertheless, if the age at which such flattening occurs decreases as the 
thickness of the lithosphere decreases, then the thinner lithosphere expected for 
Venus [Kaula and Phillips, 1981] should lead to departures from a Vx subsidence 
at a younger age than on Earth. Deviations of terrestrial seafloor depth from the Vx 
prediction that exceed the measurement errors do not occur until ~ 100 m.y. If 
comparable deviations occur for Venus at earlier ages, then this value will be an 
upper limit to the validity of the Vx relation. The distance corresponding to 100 m.y. 
age is just 100 times the apparent spreading rate in mm/y. We emphasize that this 
limit is only a guide, as there is no firm evidence to reject the notion that the Venus 
lithosphere could cool and thicken beyond this age. Moreover, even if a flattening 
of depth versus age is observed over some interval to which a Vx fit is applied, 
there is no reason to reject the thermal boundary layer model. Instead, we know 
only that a nonlinear inversion of height versus distance using the the full plate 
model or its exponential outer asymptotic solution will lead to a smaller spreading 
velocity. We conclude that, at worst, the Vx fit will lead to an overestimate of the 
spreading rate if a finite plate thickness is achieved over the length of the profile. 

For convenience, we measure elevations relative to the mean radius of Venus, 
although this may yield underestimates of the thermal boundary-layer relief if 
lowland regions with elevations of 1 km or more below the mean represent 
lithosphere that has further cooled and subsided. Any such discrepancy between 
the predicted maximum rise elevation of 1.5 km [Kaula and Phillips, 1981] and a 
greater observed relief using this datum therefore will be a minimum value. 

The second and potentially more serious source of error is the contribution to 
topography from sources other than the boundary-layer form of thermal isostasy. 
Two such causes of topography are crustal thickness variations and dynamic uplift. 
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Another possibility which could lead to large departures from a Vx relation is one in 
which thermal isostasy is the dominant mechanism of support, but the thermal 
structure is different from that expected for a simple boundary layer. Such 
situations could arise for both the vertical-tectonics and distributed-deformation 
models, wherein an overlying lid is largely decoupled from underlying flow or the 
lid moves at a velocity differing from that of the convective cell. Any of these factors 
could lead to the rejection of the hypothesis that a thermal boundary layer 
dominates the topography. 

Assuming that the spreading half-rate v is constant, the model equation relating 
distance x to elevation e for a spreading thermal boundary layer is 


e « e 0 - 8 (x/v)l/2 (4.2) 

where 8 * 350 m (m.y)' 1/2 for Earth [Parsons and Sclater, 1977] and is expected to 
be about 180 m (m.y)‘ 1/2 for Venus [Kaula and Phillips, 1981]. For Earth, we have 
performed least-squares fits to distances of 600 km and to the maximum length of 
individual profiles, solving for both eo and v. The former distance is chosen to lie 
within the known validity of the Vx subsidence for the MAR. For Venus, 600, 1500, 
and 3600 km intervals are used. The first two are chosen for comparison with 
Earth, where the 1500 km length corresponds to the mean length of the terrestrial 
profiles. The 3600 km interval spans the full topographic profiles across Aphrodite. 

Figure 4.10 shows several examples of such fits to terrestrial and venusian 
elevations. The fit to 600 km distance in Figure 4.10a, for profile 18 in the central 
north Atlantic, yields a good fit (the coefficient rof linear correlation is -0.8), and the 
rms residual topography a is modest (400 m). The spreading rate inferred from this 



140 


model is 16 ± 2 mm/y, in good agreement with the observed range of values over 
the past 80 m.y., 14-24 mm/y [Klitgord and Schouten, 1 986]. When the fit is taken 
to 2000 km, a higher a (500 m) and v (38 ± 3 mm/y) are found, because of the 
departure of subsidence from the Vx form for old seafloor and from contributions to 
topography from processes unrelated to the divergent thermal boundary layer 
(here, the Bermuda and Cape Verde rises). The residuals are uncorrelated with 
distance, except where the ends of the profile cross these rises. 

Similar fits to profile 30 in the South Atlantic are given in Figure 4.10b. The 
regression to 600 km distance yields v= 27 ± 3 mm/y with r= -0.82, and a = 300 m, 
and at 1600 km distance v= 17 ± 1 mm/y, with r= -0.88, and c = 400 m. In this 
example the apparent velocity decreases as the fit is taken farther, and the 
residuals are also uncorrelated with distance. 

Several linear regressions for e (Vx) are given in Figure 4.10c for profile 9 
across central Ovda Regio in Aphrodite Terra. The correlation is poor to 600 km 
distance (r= -0.45), but better fits (r = -0.84 to -0.86) are found for x = 1500 and 
3600 km. The large residual topography (800 m) and strong correlation of residual 
topography with distance indicates, however, that other processes contribute to the 
relief. This is not a surprising result, as several authors [e.g., Morgan and Phillips, 
1983; Banerdt, 1986; Head and Crumpler, 1987] invoke crustal thickening for Ovda 
Regio. The apparent spreading velocities are small, 1 .5 ± 0.2 and 4 ± 0.3 mm/y for 
the distance intervals 0-1500 and 0-3600 km, respectively. The corresponding 
ages of the lithosphere at the profile extrema are 1000 and 900 m.y., respectively; 
these ages so greatly exceed our 100 m.y. rule of thumb that the Vx model 
certainly may be rejected. Inferred velocities less than 1 -4 mm/y raise the question 
of just how slowly spreading can occur and still be considered plate divergence. A 
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large apparent velocity (72 ± 49 mm/y) is found for topography off of the plateau 
(1500-3600 km), and the inferred age of the lithosphere there (~20-50 m.y.) is 
within the scope of the asymptotic solution used. However, the correlation 
coefficient is very poor (-0.24) and the residuals are still correlated with x; the latter 
effect is due to the dorsae approximately perpendicular to the profile at this 
distance (see discussion of bilateral symmetry below). Such a correlation suggests 
that unmodeled processes dominate the topography, for example, crustal thickness 
anomalies [Crumpler and Head, 1988a] or local thermal basin-and-swell relief 
[Morgan and Phillips, 1983]. Moreover, it is difficult to reconcile the large 
differences in apparent spreading rate between the full profile and the distal 
segment only. 

Figure 4.10d displays several fits to profile 29, across Dali Chasma. This region 
has been suggested by Crumpler and Head [1 988a] to be more analogous to a 
typical mid-ocean ridge than the Ovda and Thetis plateaus. The data are more 
widely scattered near x * 0, yielding poor fits for both 600 and 1500 km distances. 
A reasonable fit (r = -0.66) is found over 3600 km distance, with v= 41 ±6 mm/y 
(consistent with model assumptions) and o = 500 m. As in the previous example, 
however, residual topography is strongly correlated with distance. There also 
appears to be a break in slope at ~ 1500 km distance. Taking just the interval 
beyond this break (1500-3600 km), small residuals (300 m) are found that are 
uncorrelated with distance, and r= -0.79. However, the apparent velocity is also 
small (10 ± 1 mm/y), which leads to rejection of a Vx relation. 

The rms residual topography a and coefficient of linear correlation rare given in 
Figure 4.11 for each profile as a function of regression distance interval. Figure 
4.1 la shows the results of all fits to 600 km distance. The terrestrial data are tightly 
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grouped with high (negative) rand o < 400 m; this interval near the rise crest most 
clearly shows the characteristic divergent thermal boundary topography. The 
Venus data, in contrast, are widely scattered at higher residual relief, and many 
profiles show a positive correlation of elevation with distance. In Figure 4.1 1b, the 
regressions are extended to 1500 km for Venus, and to the full length of terrestrial 
profiles, which average 1500 km. The terrestrial data are more scattered than in 
Figure 4.11a, but most values of a are less than 500 m and the magnitude of r 
exceeds 0.5. The Venus points are somewhat less scattered than in Figure 4.1 la 
but still indicate a poorer fit than for Earth. 

Figure 4.11c shows the fits to Venus only for distance intervals 0-3600 km and 
1500-3600 km. The regressions for distant topography only are nearly as poor as 
those of Figure 4.11a, so there is no evidence that only the central portions of 
Aphrodite are dominated by non-thermal boundary layer topography. The fits to 
the full 0-3600 km interval are the best in terms of r and o. Figure 4.1 Id isolates the 
values of rand o for east-central Aphrodite alone from Figures 4.1 1 a-c. Here, too, 
the only consistently good fits are those over the entire 3600 km length of the 
profiles. 

These results are summarized more compactly in Table 4.4. For the Earth, the 
majority of the variance in the topography (given by r 2 ) may be explained by the 
divergent thermal boundary model. Within 600 km of the rise crest, the fits are 
excellent, a is small (300 m), and the mean velocity and rise height match the 
observations well. Over the maximum length of the profiles, a higher rms residual 
is observed, and the velocities are overestimated due to increasing contributions 
from processes unrelated to the divergence and cooling of new lithosphere. For 
Venus, most of the fits to comparable distances are generally poor, with low rand a 
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~ 700 m. Over the maximum length of the Venus profiles (3600 km), the 'averaging 
out' of the blocky, convex profiles of parts of Aphrodite Terra leads to a higher 
correlation coefficient, but the small mean velocity (< 20 mm/y) indicates that the Vx 
relation is not valid at these distances, and so estimates of v based on the full 
solution to divergent thermal boundary layer topography would be even smaller. 
The fits to only the distant half of the profiles are very poor. For east-central 
Aphrodite alone, regressions over the full profile length show a reasonable 
correlation coefficient, modest a, a rise height close to the predicted value, and 
typical spreading rates of 30-40 mm/y. However, the goodness of fit rapidly 
degrades for other distance intervals, which suggests that the agreement may be 
fortuitous. 

Under the plate-divergence model the velocities should also show a systematic 
variation along the plate boundary due to the constraint of plate rigidity. For a 
specified angular velocity co between two plates, the velocity v as function of 
distance A from the Euler pole is 


v - cousin A (4.3) 

where R is the planetary radius. From known values of vand A, equation (4.3) may 
be solved by least squares for co; this approach depends on knowledge of 
individual plate boundaries. Figure 4.12 shows the variation of inferred spreading 
velocity with geographic location; the latter is specified simply by profile 
identification number, which increases from west to east for Venus and from north 
to south for Earth. Figure 4.12a gives the results for 0-600 km regressions for the 
MAR. The spreading rates are in the range of 10-40 mm/y, in broad agreement 
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with observations. A rotation rate of 0.32°/m.y. was recovered for the North 
American - African plate rotation by this method, which compares favorably with the 
observed value of 0.28°/m.y. [Minster and Jordan, 1 978]. Such agreement was not 
found, however, for the North American - Eurasian or South American - African 
plate pairs; in fact, the apparent velocities decrease instead of increase as A 
approaches 90°. This exercise illustrates that there is some potential to recover 
approximate relative angular velocities from topographic data alone, but the 
calculations are unreliable. 

Figure 4.12b shows the variation of v with location for Vx regressions over the 
full length of the Atlantic profiles, averaging 1500 km. The inferred spreading rates 
are larger, in the range of several tens of millimeters per year, due to the flattening 
of depth with distance for the mature thermal boundary layer and due to other 
contributions to seafloor relief. Figure 4.12c compares the apparent velocities for 
Aphrodite to 1500 km distance; in this figure zero velocity indicates a positive 
correlation of elevation with distance, for which v is undefined. Profiles with 
apparent velocities undefined or several hundred millimeters per year appear to 
separate areas with more or less regular spreading rates; inferred values are ~ 2 
mm/y in Ovda Regio, 5-10 mm/y in Thetis Regio, and 10-40 mm/y in east-central 
Aphrodite. Velocities less than 10 mm/y may not be consistent with the Vx model 
over these distances, and therefore even smaller velocities may be inferred for 
Ovda Regio, and probably for Thetis Regio as well. 

Over the full 3600 km interval in Aphrodite (Figure 4.1 2d), all apparent 
spreading rates are less than 60 mm/y and show regular variations with location. 
Again, the slowest velocities in central Ovda and Thetis Regiones are not 
consistent with the Vx form of thermal boundary layer topography. In east-central 
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Aphrodite, however, vis consistent with model assumptions and appears to 
increase systematically eastward. Such a variation might be consistent with a 
rotation pole for the this region alone located somewhere to the southwest of 
Aphrodite. The Euler pole for all CSD segments in east-central Aphrodite, 
corresponding to a time-averaged rotation pole applicable to the entire e(Vx ) fit, 
lies at 23 ± 10°N, 222 ± 30°E, in disagreement with the prediction based on 
inferred spreading rates. However, given the inability to recover accurately rotation 
rates from known terrestrial pole positions given above, this result is not conclusive. 

In summary, a divergent thermal boundary-layer component to topography 
cannot be unequivocally resolved in the spreading direction suggested by 
Crumpler and Head [1988a], and so contributions to relief other than those from a 
thermal boundary layer must generally be dominant. East-central Aphrodite alone 
is broadly consistent with spreading at v= 20-40 mm/y when averaged over the full 
profile, but the topography is much rougher than terrestrial analogs, especially near 
the rise crest. This conclusion does not exclude the plate divergence model, but 
only indicates that thermal boundary layer topography is not a discriminating 
criterion. 

TESTS FOR TOPOGRAPHIC BILATERAL SYMMETRY 

Crumpler and Head [1988a] have described both a broad, long-wavelength 
symmetry in the shape of Aphrodite and detailed matches of short-wavelength 
features on either side of the symmetry axis in a direction parallel to the CSDs. 
This mapping was based on the qualitative identification (i.e., by eye) of individual 
topographic elements in the Pioneer Venus altimetry. Crumpler and Head also 
tested quantitatively for bilateral symmetry by taking the cross-correlation between 
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the two sides of the highland, following removal of a best-fitting Vx curve. In this 
section, we apply a more detailed cross-correlation technique to the question of 
bilateral symmetry across Aphrodite Terra, and we calculate statistical confidence 
limits for these correlations. 

The cross-correlation r between two signals a and b is 

N 

(4.4) 

N N 

i=l i=l 

where N is the number of data. If the signals do not have zero mean, then the 
cross-correlation will always be nonzero. We emphasize that since ris designed to 
measure the degree of linear correlation between two variables, it is sensitive 
largely to phase differences without regard to any difference in mean amplitude. 
For example, a perfect match (r = ± 1) can be obtained between two identical 
signals differing by a multiplicative factor. In the simplest useful test of correlation, 
Pearson's r or the coefficient of linear correlation [e.g., Press etal., 1986] is given 
by subtracting the individual means of a and b before substituting into (4). For a 
topographic signal, such a test will be sensitive only to the large-amplitude, 
long-wavelength components and is appropriate for quantifying the broad bilateral 
symmetry of an elevated region. Crumpler and Head [1988a] calculated this 
cross-correlation for several profiles across Aphrodite and found values in the 
range 0.79 to 0.94. 

In order to test for bilateral symmetry of smaller elements of topography, the 
long-wavelength components must be removed. Crumpler and Head [1 988a] 
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attempted to do so by subtracting the best-fitting Vx curve from each profile; 
however, such functions lead to residual topography that is correlated with distance 
(see discussion of thermal boundary layer topography above), so it is not 
unreasonable that Crumpler and Head still found high correlation coefficients (0.65 
to 0.95). 

One approach to recovering the small-scale topography alone is through 
spectral filtering. Unfortunately, the different power spectra of Aphrodite Terra and 
the MAR preclude any simple comparison of phase coherence as a function of 
wavelength. In particular, the greater proportion of short wavelengths that must be 
phase-matched at x = 0 to form the concave and sharply peaked profile of the 
mid-ocean ridge, in contrast to the broad and often convex shape of Aphrodite, 
results in a distinctly higher cross-correlation for the Earth than for Venus. Such 
problems may be alleviated by beginning the Fourier transform at some distance 
from the rise crest, but at the expense of a loss of both spatial and spectral 
information. 

Instead, we confine our approach to the spatial domain, and we remove long- 
wavelength features by detrending the profile in a series of segments or windows. 
Each window of length L w contains N w = L w / Ax + 1 points, where Ax is the grid 
spacing. The best-fitting linear trend is removed, and the cross-correlation 
(equation 4.4), denoted by r w , is computed. These windowing operations then 
yield the variation of the cross-correlation with distance, which allows the 
recognition of individual symmetric features. We adopt a minimum window size of 
600 km, comparable to the size of the smaller features identified as symmetric 
across Aphrodite Terra by Crumpler and Head [1 988a]. An alternate size of 1200 
km is also analyzed, in order to sample better some larger features also identified 
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by these workers. 

An overall cross-correlation r p for the entire profile may be calculated from each 
point on the detrended profile. In order to make the estimate of r p more robust, and 
to ensure that any symmetric features are adequately sampled, the windows are 
overlapped by half their length. In this manner, most of the profile is doubly 
sampled, and, in effect, a second profile with detrending filters offset by L w / 2 is 
constructed. By weighting the observations inversely with the number of times they 
are sampled, the profile cross-correlation is corrected for bias introduced for the 
first half-window and last half-window of a profile which are only singly sampled 
and for the endpoints of each window which are sampled by adjacent windows. 
The origin is not included in the calculation, as its elevation is, by definition, the 
same for both sides of the profile. Finally, a global cross-correlation r g can be 
calculated in the same manner as for r p , simply by extending the summation over 
all N g data pairs in all profiles. 

The cross-correlation alone is not a complete measure of whether bilateral 
symmetry exists. Rather, a t test for rmay be formulated which gives the probability 
that the observed correlation is consistent with a sample drawn from a random 
distribution. Assuming that a and b are drawn from a joint binormal distribution, the 
test statistic for Pearson’s r is [Press etal., 1986] 

(4.5) 

and follows a t distribution with v = N - 2 degrees of freedom under the null 
hypothesis, H 0 : r = 0. The null hypothesis may be rejected at significance level a 
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in favor of the one-sided alternative hypothesis Hy; r> 0 when T> t afV . Rejection 
of this null hypothesis is the minimum quantitative requirement for the recognition 
of bilateral symmetry of specific features. Conversely, for a specified r, we may 
solve for the significance level a at which the null hypothesis may be rejected. The 
significance level represents the probability that the observed correlation is due to 
chance. A complementary expression for a is given by the confidence level (CL), 
usually expressed as a percentage: CL = 100(1 -a). For each individual window, 
then, both r w and the corresponding CL may be calculated. However, the number 
of degrees of freedom is v w = N w - 4 under linear detrending, as there are two 
degrees of freedom lost for both a and b. 

In principle, it is also possible to calculate confidence levels for both r p and r g . 
However, there is a systematic bias to the estimation procedure for these quantities 
which invalidates the corresponding null hypotheses. Although the linear 
detrending scheme is adequate to remove wavelengths much larger than the 
window size, short wavelengths which define the base or "foot" of an elevated 
region are incompletely removed by a first-order polynomial, which leads to 
correlated residuals and an enhanced correlation coefficient. 

We have calculated r g for a series of synthetic profiles in order to assess this 
effect on overall bilateral symmetry. For the Earth, a divergent-plate model 
[Parsons and Sclater, 1977] was used to generate a reference profile to a distance 
of 1500 km; we tested velocities of 10 and 20 mm/y, which approximately span the 
present spreading half-rates of the MAR. Random noise with 400 m amplitude, 
comparable to those rms residuals derived for the terrestrial profiles in the previous 
section, was added, and the global cross-correlation over 35 profiles was taken. 
This procedure was repeated ten times; the mean values for r g at v= 10 mm/y are 
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0.91, 0.36, and 0.05 for Pearson’s rand for linear detrending within 1200 km and 
600 km windows, respectively. At v= 20 mm/y, these values of r g fall to 0.88, 0.14, 
and 0.02, respectively, because the faster-spreading ridge does not flatten as 
sharply. The test was repeated for a model Aphrodite where a typical profile was 
approximated by an inverted parabola 4 km in height with a half width of 1500 km 
at the base; 38 profiles of length 3600 km were tested. The global 
cross-correlations for Pearson’s r, L w = 1200, and L w = 600 are 0.98, 0.40, and 
0.11, respectively; higher values are observed for these profiles because the 
inverted parabola has the sharpest "foot" of all three models. 

We further examined the ability to detect specific bilaterally symmetric features 
in such synthetic profiles by adding symmetric gaussian features with amplitude 1 
km (comparable in relief to those features identified by Crumpler and Head [1988a] 
and to some terrestrial oceanic plateaus) and half-width 250 km. For both 
terrestrial and venusian model topography, r g ~ 0.6-0.7 at L w = 1200 km and r g ~ 
0.3-0.4 at L w m 600 km. These values are not strongly dependent on the shape of 
the rise because the matched features dominate the cross-correlation. 

Our approach to interpreting the global cross-correlations, then, is to compare 
observations with the values expected from synthetic trials where no symmetry is 
present but which include biases similar to those present in the real data. The 
range of mean r g derived from these random models is given in Table 4.5. We 
restrict these comparisons to a qualtitative level because these models are only 
approximate; a quantitative test for the discrimination between observed and 
predicted r g is precluded since the true underlying distributions of topography are 
not well-constrained. 

Three cross-correlation tests for bilateral symmetry are applied to profile 1 8 of 

! 

! 
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the MAR in Figure 4.13. Figure 4.13a shows the two sides of the rise 
superimposed, with their individual means removed. The cross-correlation 
calculated from this approach (Pearson's r) is large (0.94) because of the 
contributions from long-wavelength topography. Figure 4.13b shows the same 
profile following linear detrending in 1200 km windows; because of the modest 
length of the profile, only two such overlapping windows are possible. The 
individual r w and CL are given in Figure 4.13c. The decay of topography 
immediately adjacent to the rise crest is adequately removed in the first window, 
and there is no significant correlation of residual topography. The second window, 
however, shows a strong signature both from the "foot" of the rise and from the 
equidistant locations of the Bermuda Rise to the west and the Cape Verde Rise to 
the east. Since both features are substantially younger then the underlying 
oceanic lithosphere [Tucholke and Mountain, 1979; Lancelot et al, 1977], the 
agreement is fortuitous. However, an argument can be made that their common 
flowline indicates a genetic relation to anomalously weak lithosphere at the rise 
crest, which facilitates the later development of hotspot volcanism [M.K. McNutt, 
personal communication, 1988]; thus the topography may be related to rise-crest 
processes. The strong correlation in the second window yields a high 
cross-correlation for the entire profile, r p = 0.73. Figure 4.13d show the same 
profile linearly detrended with L w = 600 km, giving r p = 0.26. Long-wavelength 
contributions are completely removed, and the corresponding values of r w and CL 
(Figure 4.1 3e) indicate that there is no significant bilateral symmetry at this scale. 

Profile 30 in the South Atlantic is shown for the same three calculations of the 
cross-correlation in Figure 4.14. The long-wavelength correlation (Figure 4.14a) is 
high, r p = 0.89. However, the residual cross-correlations for both the 1200 km and 
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600 km windows (Figures 4.14b and 4.14c) are negative, -0.43 and -0.42, 
respectively. In fact, Figure 4.1 4d shows that topography is negatively correlated 
for each window at L w = 600 km. 

Cross-correlations are given for profile 9 across Aphrodite Terra in Figure 4.15. 
As expected, Pearson’s ris large, 0.84 (Figure 4.15a). At L w = 1200 km (Figures 
4.15b and 4.15c), the "foot" effect of the base of Ovda Regio is clearly seen in the 
high correlation of residual topography at 1 200 km distance. The feature at x = 
3000 km is also strongly correlated and is due to the presence of symmetric ridges 
to the north and south of Ovda. These features were noted by Crumpler and Head 
[1988a] and constitute a key element in their arguments for symmetry of 
shorter-wavelength features. Because of intervening negatively correlated 
topography, however, r p = 0.10. Similar results are obtained for L w = 600 km 
(Figures 4.1 5d and 4.1 5e). 

The most strongly symmetric profile in our sample is shown in Figure 4.16. 
Profile 19 crosses central Thetis Regio and shows a consistent correlation of 
small-scale features to a distance of 2500 km for 1200 km (Figures 4.16b and 
4.16c), giving r p = 0.61. The correlation is significant ( r p = 0.49) but less impressive 
at 1^= 600 km; discarding the "foot" effect, individual matches are seen at 600 km 
and 2100 km. The nearest corresponding profile in the work of Crumpler and Head 
[1988a] shows matched features at 500, 1800, and 3000 km. 

A sample profile (number 29) across Diana Chasma is given in Figure 4.17. 
Cross-correlations for both 1200 and 600 km windows are negative, rp=-0.13 and 
-0.33, respectively. The nearest profile given by Crumpler and Head [1988a] 
shows three short-wavelength features within x = 1 200 km. Our best matches are 
at 1600 and 2700 km (Figure 4.1 7e), but both are below the 90% confidence level. 
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Profile cross-correlations are summarized by profile number in Figures 4.18 - 
4.20 for both Earth and Venus for the cases of Pearson's rand for linear detrending 
within 1200 and 600 km windows, respectively. Figure 4.18 shows generally high 
values of r p (mean ~ 0.7 for both Earth and Venus), indicating that the center of a 
broadly symmetric rise has been properly identified. For L w = 1200 km, the mean 
r p is 0.28 for Earth and 0.22 for Venus, whereas for L w = 600 km, these statistics fall 
to 0.18 and 0.13, respectively. The irregular form of Figures 4.18 - 4.20 in itself 
suggests that there is little continuity of short-wavelength features parallel to the 
rise axes. 

A better statistic of the overall symmetry is not the mean of r p , but the global 
cross-correlation r g , obtained by extending the summation of equation (4.4) over all 
profiles. These results are given in Table 4.5, along with some test values for 
synthetic profiles. Comparing the full profiles across the MAR and Aphrodite to 
these random models, the data cannot be clearly distinguished as having a higher 
correlation coefficient. We conclude that there is no strong evidence for regional 
bilateral symmetry on either planet at a scale much smaller than that of the broadly 
elevated rises. This conclusion is not changed by considering any different 
distance interval in Table 4.5. 

The absence of demonstrable regional bilateral symmetry from these tests 
conflicts not only with assertions to the contrary for Venus, but for the Earth as well. 
For the Atlantic ocean basin, Crumpler and Head [1988c] have noted several 
examples of topographic features in the Atlantic inferred to have formed at or near 
the rise crest and to have been (or to currently being) subsequently split and 
separated: Iceland, Corner Rise/Great Meteor Seamount, Ceara Rise/Sierra 
Leone Rise, and Rio Grande Rise/Walvis Ridge. We have examined 
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corresponding profiles in our data set to check if these features can be individually 
detected. Iceland itself cannot be so tested because it is presently at the rise crest; 
what few gridded data exist on the narrow Iceland ridge are negatively correlated 
(Figure 4.20). The Corner Rise and Great Meteor Seamount appear on profile 16; 
these small features are well correlated for L w = 600 km, although the Great Meteor 
Seamount dominates the pair. The correlation for the entire profile is r p = 0.30, 
which is larger than the mean value for the Atlantic (0.18) and comparable to that 
predicted for such a feature in synthetic tests (0.35). The Ceara Rise/Sierra Leone 
Rise pair do not give a good match, because the former feature is rather small and 
partially superimposed on the South American continental rise. The Rio Grande 
Rise and Walvis Ridge show a strong correlation at L w = 1200 km, and r p = 0.66; 
however, this match does not appear for the shorter window. 

We also tested five profiles on the Southeast Indian Ridge which cross the 
Kerguelen Rise and Broken Ridge. The global cross-correlations r g for these 
profiles are 0.33 for L w = 1200 km and 0.22 for L w = 600 km, and the mean values 
for r p correspond quite closely to r g . These values are lower than found for similar 
calculations for the MAR because the fast-spreading Southeast Indian Ridge (v = 
37 mm/y [Minster and Jordan, 1977]) does not show a prominent "foot" or base 
which biases correlation estimates. Synthetic profiles with a superimposed 
random component as described above show r g equal to 0.07 for L w = 1200 km 
and 0.02 for L w = 600 km. 

Therefore we find that although some individual oceanic plateaus and hotspots 
are symmetrically situated about the ridge axis, they contribute negligibly to the 
overall topography. Such a result does not preclude that similar isolated individual 
elements may exist for Aphrodite as well, but their number must be substantially 
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smaller than those features indicated by Crumpler and Head [1988a]. 

As a final test for bilateral symmetry of Aphrodite Terra in the direction 
suggested by Crumpler and Head [1988a], we examined whether a preferred 
direction of symmetry exists for profiles taken through the approximate centers of 
Ovda and Thetis. Profiles 11 and 21, respectively, pass very nearly through the 
centers of these plateaus, so a series of profiles at different azimuths was taken 
through these points. Figure 4.21a shows r p as a function of azimuth for Ovda for 
all three cross-correlation tests. The spreading azimuth suggested by Crumpler 
and Head is approximately -20°. The long-wavelength symmetry, given by 
Pearson's r, is nearly independent of azimuth, except along nearly east-west 
profiles where the presence of Thetis Regio causes a strong misfit. When 
long-wavelength components are removed using L w = 1200 km and L w = 600 km, 
the correlation as a function of azimuth varies randomly. Similar results for 
Pearson's r are found for Thetis Regio (Figure 4.21b). At both 1200 km and L w 
» 600 km, there is a slight positive correlation from azimuth -20° to +40° but no 
preferred direction. 

DISCUSSION 

We now turn to the implications of these tests for the vertical-tectonics, plate- 
divergence, and distributed-deformation models for Aphrodite Terra. We also offer 
predictions that may be tested by Magellan observations. 

The analysis of topographic bilateral symmetry has shown that if any such sym- 
metry exists, its contribution to the resolvable topography of the MAR and Aphrodite 
Terra at wavelengths of several hundred kilometers is modest at best. Since 
bilateral symmetry cannot be reliably distinguished at present, it cannot be used as 
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a test for crustal divergence. Even if the regional-scale symmetry analyzed here is 
in fact negligible, unambiguous identification of individual symmetric features (e.g., 
hotspot traces) by Magellan would nevertheless support crustal divergence, under 
either rigid-plate or distributed-deformation models. 

Our analysis of the thermal boundary layer topography of Aphrodite has shown 
that in all cases Venus scores more poorly in goodness of fit and residual 
topography than does the Earth. East-central Aphrodite alone provides the best fit 
for Venus, but under the condition that the distance of the fit is restricted to the 
longest interval. Although a thermal boundary layer component of topography 
cannot be rejected, unmodelled mechanisms contribute more than half of the 
variance in most cases. Therefore thermal boundary-layer arguments cannot be 
used to support a rigid-plate divergence model over the alternative models. 
Magellan data are unlikely to after this conclusion, as this result depends only on 
topographic length scales already resolved by Pioneer Venus. 

The burden of proof rests with the CSDs, which we have shown are not 
consistent with simple two-plate divergence. In a multiple-plate geometry, 
additional plate boundaries must be identified. The CSDs are the only identified 
major tectonic features that may be posited as such boundaries, and should be 
dominated by transform motion. In principle, if such a candidate CSD is identified 
for, say, a three-plate system, bounds on the non-transform component of motion 
(or at least its sign) across the boundary may be established by differencing the 
angular velocity vectors of the two divergent plate pairs. The angular velocity 
vectors may be constructed from Euler pole fits and topographically-inferred 
velocities. Unfortunately, uncertainties in both of these parameters preclude such 
tests. We can, however, make the following qualitative prediction: if a 



plate-bounding CSD is a pure transform, then it must lie along a small circle. In 
fact, the near-constant strike of the CSDs in Mercator projection indicates that they 
lie more nearly along rhumb lines, and so they must have either compressional or 
extensional components to their deformation. Therefore, if the plate-divergence 
hypothesis is correct, then one or more plate-bounding CSDs must exist and would 
be distinguishable from fracture zones by geological evidence for compression or 
extension across the CSD. If such evidence is not found, then the plate-divergence 
model must be rejected in favor of distributed deformation, under which CSD 
geometry is not so constrained. The vertical-tectonics model faces its strongest 
obstacle in the mere existence of the CSDs and can survive only if these 
lineaments cannot be verified. 

CONCLUSIONS 

We have considered three alternative hypotheses for Aphrodite Terra for the 
interaction of mantle convective flow with the surface: vertical tectonics, plate 
divergence, and distributed deformation. We find no convincing evidence of 
regional bilateral symmetry across Aphrodite Terra at length scales of several 
hundred kilometers, nor are we able to resolve unambiguously a thermal 
boundary-layer component to the topography of Aphrodite. Therefore neither of 
these observations may be used to support rigid-plate divergence. The existence 
of cross-strike discontinuities (CSDs) in Aphrodite implies large-scale horizontal 
motion and argues against purely vertical tectonics. The near-constant azimuths of 
CSDs are not consistent with rigid two-plate divergence and imply unusual 
apparent polar wander paths in a multiple plate system. If the plate-divergence 
hypothesis is correct, one or more CSDs must be plate boundaries and should 
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show non-transform deformation in high-resolution images. If this is not observed, 
then the CSDs are incompatible with rigid plates, and they must represent the 
traces of mantle flow upon a thin lithosphere. Accurate mapping and careful study 
of the CSDs from Magellan images will lead to the resolution of these competing 
ideas. 
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TABLE 4.1. Cross-strike discontinuity segments separating 
offset linear centers of bilateral topographic symmetry in 
Aphrodite Terra 


Latitude, 3 

°N 

Longitude, 3 

°E 

Azimuth, 

deg 

Uncertainty, 15 

deg 

-2.5 

73.0 

-22 

14 

-5.0 

77.2 

-19 

7 

-8.5 

81.5 

-23 

9 

-5.0 

94.0 

-19 

7 

-5.5 

98.7 

-19 

10 

-7.5 

107.5 

-22 

14 

-8.5 

129.2 

-26 

13 

-7.0 

137.2 

-20 

12 

-9.5 

148.2 

-26 

21 

-12.5 

156.5 

-30 

30 

-12.5 

163.5 

-38 

28 


From L.S. Crumpler (personal communication, 1987). 

a Latitude and longitude correspond to the position on the CSD 
midway between the offset centers of bilateral topographic 
symmetry. 

b Calculated from the variation in azimuth assuming that the 
uncertainty in locating individual features is ± 1° for Ovda and 
Thetis Regiones and ± 3° for east-central Aphrodite. 



TABLE 4.2. Locations of topographic profiles across Aphrodite Terra 


Number 

Latitude, 

°N 

Longitude, 

°E 

Azimuth, 

deg 

Number 

Latitude, 

°N 

Longitude, 

°E 

Azimuth, 

deg 

1 

-4.4 

62.5 

-21 

20 

-9.4 

125.0 

-21 

2 

- 3.7 

64.3 

-21 

21 

-8.8 

126.8 

-21 

3 

-2.9 

66.1 

-21 

22 

-8.6 

131.5 

-21 

4 

-2.2 

67.9 

-21 

23 

-7.9 

133.2 

-21 

5 

-1.5 

69.6 

-21 

24 

-7.2 

135.0 

-21 

6 

-7.5 

83.8 

-21 

25 

- 7.5 

140.0 

-21 

7 

- 6.7 

85.5 

-21 

26 

-6.8 

141.8 

-21 

8 

-6.0 

87.3 

-21 

27 

-6.0 

143.5 

-21 

9 

-5.2 

89.1 

-21 

28 

-16.2 

153.8 

-28 

10 

-4.5 

90.9 

-21 

29 

-15.2 

155.5 

-28 

11 

-5.0 

96.2 

-21 

30 

-12.0 

159.5 

-30 

12 

-6.2 

101.8 

-21 

31 

-16.2 

173.7 

-30 

13 

-5.5 

103.6 

-21 

32 

-15.1 

175.3 

-30 

14 

-9.4 

110.2 

-21 

33 

-14.1 

176.9 

-30 

15 

-12.8 

116.2 

-21 

34 

-13.1 

178.5 

-30 

16 

-12.1 

117.9 

-21 

35 

-14.5 

183.9 

-30 

17 

-11.4 

119.7 

-21 

36 

-13.5 

185.5 

-30 

18 

-10.8 

121.5 

-21 

37 

-12.5 

187.0 

-30 

19 

-10.1 

123.3 

-21 

38 

-11.5 

188.5 

-30 


Profiles follow great-circle paths from the specified starting point at the indicated azimuth to a 
distance of 3600 km. Profiles 4, 9, 1 1, 13, 21, 24, 25, and 29 correspond approximately to 
profiles A, G, I, J, K, L, M, and N, respectively, given by Crumpler and Head [1988a]. 
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TABLE 4.3. Locations of bathymetric profiles across the Mid- Atlantic Ridge 


Number 

Latitude, 

°N 

Longitude, 

°E 

Length, 

km 

Number 

Latitude, 

°N 

Longitude, 

°E 

Length, 

km 

1 




19 

21.3 

314.9 


2 




20 

19.2 

314.2 


3 

60.0 

330.0 


21 

17.0 

313.8 

2000 

4 

58.1 

327.3 

500 

22 

13.6 

315.4 

1400 

5 

56.1 

325.3 

900 

23 

5.2 

327.6 

1800 

6 

53.7 

324.7 

900 

24 

2.2 

329.6 

1500 

7 

50.3 

331.0 

1000 

25 

-3.8 

347.8 

1700 

8 

48.4 

332.3 

1100 

26 

-5.9 

348.8 

1700 

9 

46.5 

332.6 

900 

27 

-8.4 

346.8 

1800 

10 

44.6 

331.8 

1200 

28 

- 10.4 

346.9 

1800 

11 

42.5 

330.8 

1300 

29 

-13.2 

345.8 

1400 

12 

40.5 

330.8 

1500 

30 

-15.3 

346.3 

1600 

13 

36.3 

326.5 

2100 

31 

-17.3 

346.5 

1500 

14 

34.3 

323.1 

2000 

32 

-18.8 

347.1 

1500 

15 

32.6 

320.3 

2000 

33 

-20.7 

348.4 

1500 

16 

30.8 

318.6 

2000 

34 

-24.2 

346.9 

1500 

17 

28.5 

316.6 

2000 

35 

-27.0 

346.5 

2100 

18 

25.8 

315.3 

2000 






Profiles follow spreading flowlines (i.e., parallel to fracture zone traces) until truncated by a 
continental rise. Anomalous features include Iceland (profile 1), the Comer Rise and Great Meteor 
Seamount (profile 16), the Ceara and Sierra Leone Rises (profile 23), and the Rio Grande and 
Walvis Ridges (profile 35). 
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TABLE 4.4. Summary of least-squares tests of thermal boundary layer 
topography for Earth and Venus, e = eo + 5 (x/v) 1 / 2 


Distance Fraction Mean a, mean eo, mean v, median v, 

Interval, km with r < 0 r m km mm/y mm/y 


Earth, all profiles 


0-600 

1.0 

-0.84 

320 

4.4 

18 

16 

max. length 

1.0 

-0.80 

460 

4.1 

58 

24 

Venus, all profiles 







0 -600 

0.29 

+0.08 

660 

1.8 

39 

0 

0-1500 

0.89 

-0.47 

720 

2.8 

57 

11 

0-3600 

1.0 

-0.72 

680 

2.6 

19 

16 

1500-3600 

0.76 

-0.19 

540 

0.8 

100 

60 

Venus, East-Central Aphrodite only 





0-600 

0.28 

+0.08 

770 

0.6 

120 

0 

0-1500 

0.72 

-0.35 

590 

1.2 

20 

16 

0-3600 

1.0 

-0.66 

480 

1.3 

36 

34 

1500-3600 

0.55 

-0.23 

310 

0.7 

40 

3 




TABLE 4.5. Global Cross-Correlation Coefficients for Bilateral 
Topographic Symmetry 


Window Size, Atlantic Aphrodite Random 

km Models 


Distance Interval: full profile 


full profile 

0.70 

0.74 

0.88-0.98 

1200 

0.27 

0.17 

0.14-0.40 

600 

0.07 

0.08 

0.02-0.12 

Distance Interval: 

0-600 km 



600 

0.12 

-0.06 


Distance Interval: 

0-1200 km 



1200 

0.35 

0.20 


600 

0.12 

0.08 


Distance Interval: 

1200-3600 bn 



1200 

- 

0.14 


600 

- 

0.08 



Cumulative cross-correlation over all profiles at specified distance 
interval, taken within windows of specified length. Individual windows 
have best-fitting line removed, except where window is taken to be the 
full profile, wherein only the mean is subtracted. Random models are 
bilaterally symmetric rises with a superimposed random component 
(see text). 
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FIGURE CAPTIONS 

Figure 4.1. Topography of Aphrodite Terra, simplified from Pioneer Venus 
altimetric data [Pettengill et al., 1980; Ford ; 1986]. Mercator 
projection, 1 km contour interval. 

Figure 4.2. System of cross-strike discontinuities, or CSDs, and linear centers of 
bilateral topographic symmetry (heavy lines) in Aphrodite Terra 
mapped by Crumpler et al. [1 987] and Crumpler and Head [1 988a]. 
CSDs appear discontinuous because only the most clearly defined 
segments are included (L.S. Crumpler, personal communication, 
1 987). Topographic contours are from Figure 4.1 . 

Figure 4.3. Schematic illustrations of the three models for Aphrodite Terra tested 
in this paper; (a) vertical tectonics, (b) plate divergence, (c) 
distributed deformation. 

Figure 4.4. Best-fitting pole of rotation (Euler pole) to the azimuths of portions of 
the CSDs lying between offset linear centers of bilateral topographic 
symmetry (Table 4.1). The pole is at 20 ± 3 °N, 204 ± 17 °E (cross); 
level curves for one and two standard deviations are also shown. 
Several small circles about the Euler pole (curved lines) show the 
trajectories that CSDs would follow if Aphrodite Terra were a 
two-plate system with a rotation pole fixed with respect to the plate 
boundary. Mercator projection. 
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Figure 4.5. System of CSDs and topographic symmetry centers in Mercator 
projection about the best-fitting Euler pole. For a fixed pole 
separating two diverging rigid plates, CSDs should follow lines of 
latitude in this projection. 

Figure 4.6. Traces of topographic profiles across Aphrodite Terra analyzed in 
this paper. Each profile follows a great circle approximately 
perpendicular to the linear center of bilateral topographic symmetry. 
Profile spacing is approximately 200 km, and centers are taken no 
closer than about 150 km to the nearest mapped CSD. Profiles cross 
distant CSD segments in only a few cases. 

Figure 4.7. Traces of bathymetric profiles across the Mid-Atlantic Ridge (MAR), 
for comparison with Aphrodite Terra. Traces follow spreading 
fiowlines [Klitgord and Schouten, 1986; Cande etal., 1988] in most 
cases, except for short profiles on the Reykjanes Ridge and near the 
equator, which are taken along small circles about the best-fitting 
two-plate rotation pole [Minster and Jordan, 1 978]. Profile centers 
are selected in a similar fashion to those for Aphrodite Terra, spaced 
approximately 200 km apart and no closer than about 150 km to 
large-offset (> 1°) transform faults. Profiles are truncated where they 
intersect continental rises or oceanic plateaus with continental 
basement. 
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Figure 4.8. Topographic profiles of Aphrodite Terra constructed from individual 
Pioneer Venus measurements. Figure is divided into Ovda Regio 
(60-110° E), Thetis Regio (110-150° E), and east-central Aphrodite 
(150-190° E). Within each region, the longitude of the profile center 
increases upward on the figure. Profiles are offset by 3 km in 
elevation for clarity. 

Figure 4.9. Bathymetric profiles across the Mid-Atlantic Ridge, constructed from 
DBD5 digital bathymetry degraded to Pioneer Venus resolution (see 
text). Figure is divided by plate pair (North American - Eurasian, 
North American - African, South American - African) for convenience. 
For each plate pair, profiles are arranged from north to south going 
down the figure, and are offset by 3 km in depth for clarity. 

Figure 4.10. Least-squares fits of elevation e versus the square root of distance x 
from topograpic profiles across Aphrodite Terra and the Mid-Atlantic 
Ridge. Statistical parameters and apparent rates of divergence are 
given in the text for each fit. 

(a) MAR profile 18 (north of the Kane Fracture Zone); fits to 600 and 
2000 km distances. 

(b) MAR profile 30 (north of the St. Helena Fracture Zone); fits to 600 
and 1 600 km distances. 

(c) Aphrodite Terra profile 9 (Ovda Regio); fits to 600, 1500, 3600 km 
distances and to the interval 1500-3600 km. 

(d) Aphrodite Terra profile 29 (east-central region); fits to 600, 1500, 
3600 km distances and to the interval 1500-3600 km. 
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Figure 4.11. Summary of rms error o versus linear correlation coefficient r fore 
versus Vx regressions for all profiles on Earth (circles) and Venus 
(triangles). 

(a) Regressions to 600 km distance. 

(b) Regressions to 1500 km distance for Venus and full profile length 
(averaging 1500 km) for Earth. 

(c) Regressions for Venus only to 3600 km distance (triangles) and 
for the distance interval 1500-3600 km (squares). 

(d) Regressions for east-central Aphrodite alone to distances of 600 
km (squares), 1500 km (circles), 3600 km (triangles), and for the 
interval 1 500-3600 km (crosses). 

Figure 4.12. Variation of apparent half spreading rate vwith geographic location, 
given by profile number. 

(a) MAR, 600 km fits. The boundary between the Eurasian and 
African plates occurs between profiles 12 and 13, and that between 
the North American and South American plates is between profiles 
21 and 22. 

(b) MAR, full length of individual profiles (mean 1500 km). 

(c) Aphrodite, 1500 km fits. Zero velocity indicates a positive 
correlation of elevation with distance; for such profiles a 
measureable cooling boundary layer component may be rejected. 

(d) Aphrodite, 3600 km profiles. 
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Figure 4.13. Tests of bilateral symmetry for MAR profile 18. 

(a) West (solid) and east (dashed) sides of profile 18, with means 
subtracted. 

(b) Same profile with windowed detrending operations applied to 
remove long-wavelength components of topography. Window length 
in this example is 1200 km, and the best-fitting line has been 
removed. Overlapping windows offset by 2 km in depth for clarity. 

(c) Correlation coefficient (solid line) as a function of distance for the 
windowed, detrended profiles in (b). Dashed lines show confidence 
levels at which the null hypothesis for each window may be rejected. 

(d) Atlantic profile 18 linearly detrended within 600-km windows. 
There are now five individual, overlapping window segments within 
the profile. 

(e) Correlation coefficient as a function of distance lor (d). 

Figure 4.14. Tests of bilateral symmetry for MAR profile 30. 

(a) Full profile to 1600 km, with means subtracted from each side. 

(b) Profile after linear detrending in 1200 km windows; since the 
profile length is1600 km, only a single window is used. 

(c) Four overlapping 600-km windows for the same profile. 

(d) Correlation coefficients and confidence levels for (c). 

Figure 4.15. Tests of bilateral symmetry for Aphrodite Terra profile 9, Ovda Regio. 

(a) Full profile (3600 km), with means subtracted. 

(b) Profile linearly detrended in 1200-km windows. 

(c) Correlation coefficients and confidence levels for (b). 
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(d) Profile linearly detrended within 600-km windows 

(e) Correlation coefficients and confidence levels for individual 
windows for(d). 

Figure 4.16. Tests of bilateral symmetry for Aprodite Terra profile 19 (Thetis 
Regio), the profile observed to be the most highly symmetric at short 
wavelengths. 

(a) Full profile (3600 km), with means subtracted. 

(b) Profile linearly detrended in 1200-km windows. 

(c) Correlation coefficients and confidence levels for (b). 

(d) Profile linearly detrended within 600-km windows 

(e) Correlation coefficients and confidence levels for individual 
windows for(d). 

Figure 4.17. Tests of bilateral symmetry for Aphrodite profile 29 (East-Central 
Region, across Dali Chasma). 

(a) Full profile (3600 km), with means subtracted. 

(b) Profile linearly detrended in 1200-km windows. 

(c) Correlation coefficients and confidence levels for (b) 

(d) Profile linearly detrended within 600-km windows. 

(e) Correlation coefficients and confidence levels for individual 
windows for(d). 


Figure 4.18. Summary of profile correlation coefficients, given by profile number, 
for simple cross-correlations of full profiles. Strong correlation is 
evident in all cases for this method because of the broad central 
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elevated region, (a) Earth, (b) Venus 

Figure 4.19. Summary of profile correlation coefficients under linear detrending 
within 1200-km windows. On both planets, systematic contributions 
to cross-correlations arise from correlated residuals after linearly 
detrending across the base of elevated region, (a) Earth, (b) Venus. 

Figure 4.20. Summary of profile correlation coefficients under linear detrending 
within 600-km windows. There is no consistent evidence for 
symmetric features at this scale on either planet, (a) Earth, (b) 
Venus. 


Figure 4.21. Profile correlation coefficients as a function of azimuth from the 
approximate centers of (a) Ovda Regio and (b) Thetis Regio. Solid 
line, mean topography subtracted; long dash, linear detrending 
within 1200-km windows; short dash, linear detrending within 600 km 
windows. 
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CHAPTER 5: 

PENECONTEMPORANEOUS METAMORPHISM, FRAGMENTATION, AND 
REASSEMBLY OF ORDINARY CHONDRITE PARENT BODIES 
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INTRODUCTION 

Studies of the thermal evolution of ordinary chondrites are founded on obser- 
vations of two basic parameters: (1) the distribution of peak metamorphic 
temperatures for the various petrologic types, as inferred from mineralogy, and (2) 
cooling rates as measured by metallographic, fission track, or radiogenic gas 
retention methods. The correlation between these parameters is not clear, mainly 
due to disagreement between the various cooling rate methods. Simple models for 
internally (i.e., radioactive ly) heated chondrite parent bodies should produce an 
inverse correlation, that is, slow cooling in the deep interior for higher petrologic 
types, and more rapid cooling near the periphery for lower petrologic types. 
However, Scott and Rajan [1981] noted that there appears to be no correlation 
between the metallographic data and petrologic type, so they devised a 
"metamorphosed-planetesimal" model to resolve this discrepancy, in which heating 
occurs in small objects that are subsequently randomly accreted into larger parent 
bodies. In this paper, I review the thermal and accretional requirements of the 
metamorphosed planetesimal model and suggest an alternative, consistent with 
the metallographic cooling rate constraints, in which ordinary chondrite parent 
bodies are collisionally fragmented and then rapidly reassembled before 
metamorphic heat has been dissipated. 

THERMAL HISTORY OF ORDINARY CHONDRITES 

Ordinary chondrites have been heated to temperatures varying from ~400°C for 
type 3 to ~950°C for extreme type 6 [Dodd, 1981]. The classification of ordinary 
chondrites by chemical group (H, L, LL) is not of direct importance here except to 
constrain the different groups to separate parent bodies or accretion zones with 
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little mixing between them. The volume abundances of the various petrologic types 
[Wasson, 1974], although different for each group, show that the more intensely 
metamorphosed type 5 and 6 materials are the most common. Reproduction of 
petrologic type distribution, and hence distribution of maximum metamorphic 
temperatures, is the first constraint that thermal models of chondrite parent bodies 
must satisfy. 

Cooling rates provide a second constraint on thermal models. Cooling rates in 
ordinary chondrites have been measured directly by metallographic and 
fission-track methods and have been inferred indirectly from isotopic ages. The 
metallographic method [Wood, 1967], using Ni-diffusion in metal phases, yields 
cooling rates of 1-10°C/m.y. as the material cooled through ~500°C (Ni-diffusion 
becomes sluggish at lower temperatures). An important result of these 
measurements, pointed out by Scott and Rajan [1981], is that there appears to be 
no correlation between cooling rate and petrologic type in the 27 measurements 
summarized by Wood [1979]. Although simple models, such as the "onion-shell" 
model discussed below, should produce some correlation, it is not clear that the 
metallographic measurements should always reflect this for lower petrologic types. 
That is, such measurements of cooling rates near 500°C for type 3 chondrites may 
reflect initial slow cooling just beyond the peak temperature of these meteorites, 
while similar measurements of type 4-6 chondrites represent relatively rapid 
cooling well below peak temperatures of these materials. Hence the issue of 
correlation betweeen metallographic cooling rate and petrologic type may not be 
meaningful if the cooling rate measurements reflect different parts of the 
metamorphic history of, say, type 6 and type 3 chondrites. However, in spite of the 
fact that I will exclude type 3 chondrites from this analysis, the lack of correlation of 
petrologic type with metallographic cooling rate among types 4-6 still persists. 
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Furthermore, other cooling rate data are in disagreement with the metal- 
lographic method. Pellas and Storzer[ 1979] obtained H chondrite cooling rates 
from 244 Pii fission tracks that are inversely correlated with petrologic type 
(~3-7°C/m.y. for type 6, 7-17°C/m.y. for type 5, and >10-44°C/m.y. for type 4). 
Similar measurements on L chondrites [Pellas, 1982] are more perplexing. Some 
have cooling rates inversely correlated with petrologic type, as in the H chondrites, 
except slower than corresponding metallographic rates. Others have extremely 
rapid cooling rates ("almost instantaneous" within the resolution of the method) for 
L5-6 material and slower rates for L4-5 samples, thus showing a direct correlation 
between cooling rate and petrologic type, the reverse of the H case. 

Turner et al. [1978], using 40 Ar- 39 Ar methods, deduced that ordinary chondrites 
must have passed through the argon closure temperature interval (360-1 00°C) 
between 50 and 100 m.y. after formation. Wood [1979] has noted that metallo- 
graphic cooling rates are slower than Ar retention ages by a factor of about six not 
only for ordinary chondrites, but also for irons, pallasites, and mesosiderites. Willis 
and Goldstein [1981] have revised the original ordinary chondrite metallographic 
cooling rates upward by a factor of two. 

In summary, the variations in measured cooling rates and correlation with 
petrologic type between different methods preclude a simple unified picture. 
Diverse techniques yield direct correlations (some L group fission track), inverse 
correlation (fission track for H and some L group), or no correlation (metallograph- 
ic). For the rest of this discussion, I will assume that any error in the metallograph- 
ic method is a constant factor for type 4-6 chondrites and therefore that the lack of 
correlation between metallographic cooling rate and petrologic type is significant. 
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ONION-SHELL MODEL 

Conventional ("onion-shell") thermal models of ordinary chondrite parent 
bodies [e.g., Herndon and Herndon, 1977; Minster and Allegre, 1979; Miyamoto et 
al., 1981] envision rocky asteroids of ~100-km radius that are accreted cold in ~10 5 
yr and are subsequently heated by decay of 26 AI (half-life 7.2 x 10 5 yr). Other heat 
sources have been proposed, such as other short-lived radionuclides, superheavy 
elements, electrical induction, and solar, tidal, and impact heating. Sonett and 
Reynolds [1979] have reviewed the problem of primordial heat sources. Since the 
discovery of relict 26 Mg in Allende refractory inclusions by Lee et al. [1976], 26 AI has 
been the favorite such heat source. Evidence for 26 AI has been strengthened 
recently by the detection of relict 26 Mg in an anorthite-bearing chondrule in the 
Semarkona chondrite [Hutcheon et al., 1988]. For consistency with previous 
models, 26 AI is the assumed heat source for the rest of this discussion. The princi- 
pal conclusions of this work are not affected by the specific choice of heating mech- 
anism, as long as the time scales are similar (e.g., solar electromagnetic heating). 

The onion-shell model requires rapid accretion of bodies several hundred 
kilometers in diameter prior to significant 26 AI decay (neglecting the formation 
interval of 26 AI relative to the formation of the solar system). This time scale (~10 5 
yr) is consistent with numerical simulations of the collisional and gravitational 
interactions of early solar system bodies of Greenberg et al. [1978]. On the other 
hand, the theory of gas-free accumulation [Safronov, 1979; Wetherill, 1980] 
predicts accretion time scales of >10 7 yr for such bodies in the vicinity of the 
asteroid belt. This latter time scale is a necessary condition for the 
metamorphosed-planetesimal model, discussed below. Since the problem of 
accretion remains moot, assumptions of accretion time scale in the context of 
chondrite thermal models must be viewed as ad hoc. 
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The onion-shell models produce type 6 material in their cores, surrounded by 
successive zones of types 5, 4, and finally type 3 material near the surface (Figure 
5.1). The observed distribution of petrologic types among H chondrites can be 
produced using a parent body with a radius of 85 km, thermal diffusivity k of 5x1 0* 3 
cm 2 s* 1 , and initial 26 AI/2 7 AI ratio of 5x1 0* 6 [Miyamoto et a/., 1981]. Cooling rates are 
in agreement with fission track measurements by Pellas and Storzer [1979] and are 
inversely correlated with petrologic type. An inverse correlation is expected for 
such a model, as regions near the core cool slowly and material nearer to the 
periphery cools more rapidly. 

In the context of metallographic cooling rates, however, the onion-shell model 
fails distressingly, as these data indicate no correlation with petrologic type. Hence 
some two-stage process apparently is required to fit the metallographic data, first 
metamorphosing the ordinary chondrites, and then transporting them to random 
cooling locations prior to slow cooling through 500°C. Two such models will be 
described and compared. 

METAMORPHOSED-PLANETESIMAL MODEL 

A "metamorphosed-planetesimal" model has been proposed by Scott and 
Rajan [1981] and Taylor et at. [1982]. They suggest that if accretion of 100-km- 
radius bodies takes £10 7 yr and if 26 AI decay is the dominant heat source, then 
maximum temperatures are reached (i.e., petrologic types are established) in small 
planetesimals with r<, 10 km. Cooling rates in these bodies are controlled by their 
burial depth as they are subsequently randomly accreted into larger r~100 km 
parent bodies; therefore cooling rates appear uncorrelated with petrologic types 
(Figure 5.2). 

This model rests on several assumptions that must be explored more 
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thoroughly. Consider the first stage of the process, the heating of r £ 10 km 
planetesimals. In order to heat these objects sufficiently, they must initially have 
either a greater proportion of radionuclides or a much poorer thermal diffusivity 
than r ~100 km bodies, or both. Wood [1979] describes such a model, wherein 
particulate matter (k=10‘ 5 cm 2 s -1 , 26 AI/2 7 AI = 6x1 0- 5 ) accretes at exponentially 
decreasing rates. Upon reaching a critical temperature, the "dust" grains 
irreversibly sinter to rocky, chondritic material (K=5x10‘ 3 cm 2 s- 1 ). Wood's model 
produces melting temperatures on bodies that grow to 10-20 km radius. 

Yomogida and Matsui [1984] have developed an improved sintering model, but 
their results differ from Wood's largely because they do not consider accretion and 
use only long-lived radionuclides. An interesting feature of their model is that peak 
temperatures are nearly uniform in the interior, more so than even the onion-shell 
model, which produces such large proportions of petrologic types 5 and 6. In the 
sintering model, large gradients exist only at the unsintered, powdery surface layer. 
This implies that planetesimals of different sizes are each composed of nearly a 
single petrologic type, in agreement with the assumptions of the metamorphosed- 
planetesimal model. 

Hence a major difference between the metamorphosed-planetesimal and the 
onion-shell parent body is the consolidation state of the object in which 
metamorphism occurs. The onion-shell model is based on the assumption that 
solid bodies exist before significant endogenic heating begins (10 5 -10 6 yr). 
Carbonaceous chondrites, for example, have not been heated above ~300°C, yet 
many are highly compacted [Cain and McSween, 1984]. Since the physical 
properties of primitive solar system material remain unknown, neither model can be 
excluded. Material properties are also important in collisional studies as discussed 
below. 
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The second step in the metamorphosed-planetesimal model is the accretion of r 
<10 km bodies into r~100 km objects on a timescale £ 10 7 yr. An implicit 
assumption here is that the planetesimals accrete fairly intact (i.e., like a bunch of 
grapes) in order to retain sufficient heat so that cooling rates at 500°C can only be 
achieved under one of the following circumstances: (1) fragments are rapidly 
covered by other impacting debris, (2) fragments are emplaced by impact at a 
depth of -4 km under regolith or -40 km under rock, or (3) planetesimals are 
shattered but not dispersed by impact, in a manner analogous to a snowball 
striking a wall (R. Greenberg, personal communication, 1984). The required long 
accretion times of 100 km bodies implies that fragments will not be covered by 
debris from other incoming planetesimals quickly enough to prevent rapid cooling. 
On the other hand, the maximum depth of emplacement can be estimated by 
considering the impact of a 10-km radius planetesimal at 100 m/s (the escape 
velocity of a 3 g/cm 3 , 77-km radius body). A crater diameter of -13 km is predicted 
in a crystalline target [Dence et al., 1977]. A typical diameter: depth ratio of 5 for 
craters of this size [Gault et al., 1975] gives a crater depth of 2.6 km. Earlier 
accreting planetesimals will form smaller and shallower craters because of lower 
impact velocities. Thus even though this calculation is strictly valid only for 
hypervelocity impacts onto surfaces of planet-size objects, it is apparent that 
isolated fragments of a shattered metamorphosed planetesimal can cool slowly 
only if they are completely buried at maximum excavation depths. Since it seems 
unlikely that much material was emplaced this way, I conclude that the metamor- 
phosed-planetesimal model is valid only if the r< 10 km bodies can accrete more 
or less intact into r -100 km chondrite parent bodies or if the fragmentation of 
accreting planetesimals is strongly nondispersive. Models of the collisional 
evolution of objects in the solar system can illuminate this problem. I have tested 
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the metamorphosed-planetesimal accretion scheme against such a model and 
have used it to search for other mechanisms to explain apparently uncorrelated 
metamorphic grades and metallographic cooling rates in ordinary chondrites. 

COLLISION DYNAMICS 

A significant result of the work of Greenberg et al. [1978] applicable to the 
present problem was the development of a comprehensive collision model. Impact 
velocity is given by Vj 2 =v 0 2 + v e 2 , where v 0 is the encounter velocity, and v e is the 
mutual escape velocity of the two bodies. Outcomes of collisions are divided into 
four categories: (1) rebound intact, (2) cratering of both bodies and rebound, (3) 
smaller body shattered, larger cratered, and (4) both bodies shattered. The result 
for a particular collision depends on the impacting bodies' resistance to both 
simple cratering and catastrophic disruption. These results are accompanied by 
the accretion of none, part, or all of the impacting material, depending on the 
partitioning of kinetic energy between the two impactors and their ejecta. Data for 
materials ranging from loosely bound regolith to basalt were included. 

The collision algorithm of Greenberg etal. [1978] may be used to investigate the 
effects of accretion mechanics on metamorphic models of ordinary chondrites, 
employing parameters given by these workers for solid rock and for an 
intermediate material between rock and regolith (perhaps corresponding to coher- 
ent dirt clods). These parameters were chosen as representative endmembers for 
material properties of early solar system planetesimals or meteorite parent bodies. 
Greenberg et al. [1 978] give full details of the method. 

Two characteristic encounter velocities were studied. Safronov [1972] analyt- 
ically calculated that equilibrium encounter velocities are comparable to the 
escape velocities of the largest bodies in the population. Accordingly, the upper 



218 


velocity used was 1 00 m/s (velocities of several hundred m/s may be applicable but 
do not affect the conclusions). On the other hand, Greenberg et al. [1978] found 
that encounter velociites are comparable to the escape velocities of the original 
bodies in the population, i.e., a few m/s. Hence a lower velocity of 5 m/s was 
adopted. 

The lower velocity collision results are shown in Figure 5.3. The most important 
feature to note is that, except for a small field of nearly equal-sized rocky 
planetoids, collisions involving a body £ 50 km in diameter cannot result in 
accretion intact or with relatively minor cratering. In terms of the metamorphosed- 
planetesimal model, this means that impacting 10 km planetesimals will be 
shattered when the parent body has grown to only 20-30-km radius. Heat losses 
for dispersed fragments will be rapid compared with the required long accretion 
timescale, so the metamorphosed-planetesimal model breaks down in this case, 
unless the fragments remain closely aggregated. 

At higher encounter velocities (Figure 5.4), there are no fields at all that allow 
"intact” accretion. All outcomes result in the shattering of one or both bodies, and 
objects £ 50 km do not have sufficient gravitational binding energy to reaccrete; 
instead the fragments escape. Increased encounter velocities result in progres- 
sively more disruption. Therefore the metamorphosed-planetesimal model is 
probably most applicable to strong, rocky bodies at low encounter velocities, so 
that accretion proceeds without shattering. Since high impact strength is needed in 
this case, metamophosed- planetesimals must be more or less completely sintered 
to rock before they are accreted. Even then, r -100 km objects must be accreted in 
a stepwise fashion, first assembling several r<10 km planetesimals and then 
accreting several of these intermediate aggregates (so as to pass from field I to 
field II in Figure 5.3b). These constraints may be relaxed if nondispersive 
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shattering of impacting planetesimals occurs. 

FRAGMENTATION-AND-REASSEMBLY MODEL 

Is there an alternative accretion scheme to explain uncorrelated petrologic 
types and cooling rates? Figures 5.3 and 5.4 show sizable areas (field IV) where 
the kinetic energy delivered per unit volume is sufficient to catastrophically 
fragment both bodies, but the remains will reaccrete into a gravitationally bound 
"rubble pile" [Davis et al., 1979]. At the relatively low encounter velocities 
described here, this result occurs for comparable size bodies (larger: smaller 
diameter ratio < 5), with the larger object's diameter £100 km. Hartmann [1979 a,b] 
has investigated collisions of comparable size bodies for a variety of initial 
conditions. He concludes that impacts between comparable size bodies may 
account for several peculiar asteroid observations by producing unfractured or 
partially brecciated contact binaries, totally brecciated spheroids, or coorbiting 
groups. For example, Trojan asteroid 624 Hektor may be a highly fractured contact 
binary. 

A series of collisions between comparable-size objects could rescue the meta- 
morphosed-planetesimal model, as described at the end of the previous section. 
Such collisions could also account for the inferred ordinary chondrite thermal 
histories in the following manner. 

First, rapid accretion causes peak temperatures due to 26 AI heating to occur in 
objects of ~100 km radius, as in the onion-shell model. The relative proportions of 
the various petrologic types are established at this stage (Figure 5.5a). 

Second, maximum metamorphic temperatures are maintained for several tens 
of million years in bodies of this size. Collision with another object on this time 
scale completely fragments both bodies (Figure 5.5b). Collisions capable of 
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shattering a hypothetical onion-shell parent body are comparatively rare in the 
present solar system. Using the collision model of Greenberg etal. [1978] and the 
impactor fluxes given in Greenberg and Chapman [1983], and assuming a mean 
encounter velocity of 5 km/s, I calculate that such collisions should occur 
approximately on a time scale varying from 10 8 to 10 10 yr, depending on target 
material properties. In contrast, the earliest planetesimal evolution may have been 
characterized by extremely common, mutually fragmenting encounters. A 
characteristic fragmentation encounter time for r ~100 km objects of only 10 2 to 10 4 
yr is predicted from the numerically simulated planetesimal distributions produced 
at ~2x1 0 4 yr by the Greenberg et al. [1978] model. Hence a fragmentation 
timescale of ~10 7 -10 8 yr, as required by the characteristic cooling time of the 
onion-shell model, is not inconsistent with plausible planetesimal size and speed 
distributions intermediate between the earliest and the present-day populations. 

Third, the fragments reassemble into a megaregolith rubble pile while they are 
still hot (Fig. 5c). Factors affecting reassembly time are complex (such as size and 
speed distribution of fragments and partitioning of translational versus rotational 
energy), but an approximate estimate can be made assuming that gravitationally 
bound fragments reaccrete promptly and do not form coorbiting swarms. Under 
this assumption, a characterstic reaccretion time is given by the ballistic flight time 
of ejecta. In the Greenberg et al. [1978] model, the velocity of fragments is 
computed such that all have the same speed and their kinetic energy is 50% of the 
impact energy, an upper limit for their experiments with basalt. For two 3 g/cm 3 - 
100-km radius asteroids with an encounter velocity of 100 m/s, the characteristic 
fragment velocity is 135 m/s (0.83V e ), and the ballistic flight time is about 3.5 hr. 
Therefore fragments greater than a few meters in diameter could retain sufficient 
heat so that cooling rates at 500°C are controlled by burial depth in the reaccreted 
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parent body. For example, a 10-m diameter fragment would require around 45 
days to cool from 800°C to 500°C at its center. 

The likelihood of a significant fraction of the shattered chondrite parent body 
occurring in such sizeable fragments can be estimated from a mass distribution 
model. Greenberg et al. [1978] used a conventional power-law mass distribution of 
the form N = Cm - b where N is the number of fragments with mass > m, and C and b 
are constants. Hartmann [1969] found that b varies from about 0.67 when 
fragmentation energy is minimal to near unity in high-energy collisions or where 
extensive regrinding takes place. For a fragmented 100-km radius parent object, 
and for values of b up to ~0.8, over 99% of the mass is found in fragments >1 0 m in 
diameter. Thus even for more energetic collisions (although they might preclude 
reassembly) or encounters in which secondary collisions cause further 
comminution, much of the mass is distributed in fragments capable of retaining 
metamorphic heat for many days, perhaps months. If these fragments are 
reaccreted before much of this time has elapsed, then cooling rates will appear 
uncorrelated with petrologic type, assuming that the fragments are thoroughly 
mixed by the collision. 

DISCUSSION 

The fragmentation-and-reassembly model of chondrite parent bodies requires 
two of the same assumptions as does the metamorphosed-planetesimal model. 
First, accretion must preserve overall homogeneity among H, L, and LL chondrites; 
that is, there is little interaction between these objects (or groups of objects). 
Second, the accreting objects must have comparable thermal states, i.e., heated to 
ordinary chondrite metamorphic temperatures ~400-950°C. If not, rapid heat 
losses to cold fragments would cause cooling rates much larger than the usual tens 
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of degrees per million years or so (however, recall that some L chondrites have 
very fast cooling rates). 

Other workers have suggested that meteorite parent bodies may have been 
collisionally fragmented and then gravitationally reassembled one or more times 
over solar system history [Hartmann, 1979 b; Taylor et al., 1982]. Greenberg and 
Chapman 1983] calculated that asteroids in the diameter range 30-300 km have 
undergone this process several to many times (the number of such events 
decreases with larger objects). However, none consider the possibility of such 
collisions occurring early, while parent bodies were still at metamorphic 
temperatures. 

Rubin etal. [1983] found that some fragmental breccia (predominantly L and LL) 
had incorporated melt rock clasts prior to metamorphism, while others (mostly H) 
had acquired melt rock clasts after metamorphism. Scott and Rajan [1 981] studied 
metallographic cooling rates in xenolithic clasts in four ordinary chondrite breccias. 
Cooling rates in clasts from two of the meteorites varied widely, whereas clasts in 
the other two had coherent cooling rates. The former case was interpreted as 
compaction of host and clasts after cooling through 500°C, and the latter as 
incorporation of clasts before or during metamorphism. Scott and Rajan [1 981] and 
Taylor et al. [1982] invoke fragmentation and reassembly to explain association of 
breccia clasts with a wide range of cooling rates (the resulting wide range of burial 
depths could probably not be sampled without shattering the parent body.) 
However, they propose metamorphosed-planetesimals to account for the lack of 
correlation between petrologic type and cooling rate and for penecontempora- 
neous compaction and metamorphism of clasts in other meteorites. Repeated 
fragmentation and reassembly of the ordinary chondrite parent bodies is a single 
mechanism that can explain both uncorrelated cooling rates and petrologic types 
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and incorporation of melt-rock and xenolithic clasts before, during, or after 
metamorphism. 

CONCLUSIONS 

The metamorphosed-planetesimal model has been proposed previously to 
explain the apparent lack of correlation between petrologic type and metallo- 
graphic cooling rate by suggesting that if accretion of r~100 km bodies takes >10 7 
yr, then peak temperatures due to 26 Al heating are reached in objects of r£lO km. 
Cooling rates at 500°C are then fixed by burial depth of these planetesimals in the 
growing parent body. The model requires that initial materials have much lower 
thermal diffusivities and/or greater proportions of radionuclides than assumed for 
conventional onion-shell parent bodies. In addition, the metamorphosed- 
planetesimal model demands relatively intact accretion or nondispersive 
fragmentation of r<10 km bodies so that cooling rates are controlled by burial 
depth. If these requirements are not met, fragments will cool rapidly near the 
surface instead of slowly at depth. 

Alternatively, onion-shell parent bodies could be completely fragmented and 
then gravitationally reassembled by low-velocity collision of comparabe-size 
objects. If such collisions occurred during the few tens of m.y. of peak metamorphic 
temperatures, then cooling rates at 500°C would be determined by burial depth in 
the reaccreted parent body, because reassembly times may be quite short. 

All models described here are chiefly dependent on the thermal properties and 
radionuclide proportions of primitive solar system material, the style and timescale 
of accretion, and the "true" nature of the relationship between cooling rate and 
petrologic type. No model can be eliminated until the answers to these questions 
are better resolved. 
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FIGURE CAPTIONS 

Figure 5.1: Sample "onion-shell" rocky parent body for H chondrites showing 

successive zones of petrologic types 3-6, after Miyamoto et al. [1981]. 
Peak temperatures are fairly uniform in the deep interior and drop 
smoothly near periphery. Cooling rates at 500°C are predicted to be 
inversely correlated with petrologic type. 

Figure 5.2: "Metamorphosed-planetesimal" model, after Scott and Rajan [1981]. 

Chondrite metamorphism occurs in small, loosely consolidated 
planet- esimals, which then accrete into larger parent bodies. 
Cooling rates at 500°C are controlled by burial depth and appear 
uncorrelated with petrologic type. 

Figure 5.3: Results of low velocity two-body collisions using model of Greenberg 

et al. [1978]: l-bodies accrete intact; ll-bodies crater each other and 
accrete; Ill-smaller body shattered, larger cratered, fragments accrete; 
IV both bodies shattered, fragments accrete; V-as result III, fragments 
escape; Vl-as result IV, fragments escape. Brackets on ordinate and 
absicissa give approximate sizes of metamorphosed-planetesimals 
and parent bodies into which they accrete, respectively. Metamor- 
phosed-planetesimal is valid under outcome I or II; applicability under 
result III is doubtful (see text). Result IV shows fragmentation and 
reassembly, providing an alternative model. 
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Figure 5.4: As Figure 5.3, but with higher encounter velocity. Smaller bodies 

have insufficient gravitational binding energy to prevent disruption 
following collision, and escape fields dominate as velocities approach 
1 km/s. 

Figure 5.5: Fragmentation and reassembly model. Chondrite metamorphism 

occurs in 100-200-km diameter bodies that are fragmented by 
collision with a comparabe-size object and then gravitationally 
reassembled within a few hours-days. Cooling rates for fragments 
greater than a few meters in diameter are controlled by burial depth in 
the reaccreted parent asteroid. 
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(b) Rocky material 
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CHAPTER 6: 


WATER AND THE THERMAL EVOLUTION 
OF CARBONACEOUS CHONDRITE PARENT BODIES 
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INTRODUCTION 

Carbonaceous chondrites, though among the most chemically primitive mate- 
rials in the solar system, have been nonetheless modified within their parent 
bodies. Aqueous alteration, in particular, has played an important role in the 
petrogenesis of CM and Cl chondrites [Kerridge and Bunch, 1979; McSween, 
1979]. In the first attempt to understand the link between thermal evolution and 
aqueous alteration, DuFresne and Anders [1 962] proposed a qualitative model for 
a carbonaceous chondrite parent body in which internal heating released water 
from pre-existing hydrous minerals until a liquid zone extended to just below the 
surface. Aqueous alteration was thought to have occurred there over a period of a 
few thousand years, with the retention of free H 2 0 on the parent object limited by 
subsolar sublimation of ice. More recent research on carbonaceous chondrites 
has not addressed the thermal driving mechanism for alteration, but has focused 
on the complex mineralogy of altered matrix and vein materials [Richardson, 1978; 
Bunch and Chang, 1980; Barber, 1981 ; McKinnon and Zolensky, 1984; Tomeoka 
and Buseck, 1985, 1988; McSween, 1987]. Oxygen isotope systematics [Clayton 
and Mayeda, 1984] have nevertheless provided constraints on temperature and 
water volume present during the alteration process. 

In this paper, we develop the first quantitative models for the thermal evolution 
of carbonaceous chondrite parent objects in which the thermodynamic and 
transport properties of water are shown to exert -a strong influence. In contrast to 
the approach taken by DuFresne and Anders [1962], we assume that these parent 
bodies were initially mixtures of ice and anhydrous rock. We argue that there are 
two likely alternatives for the parent body aqueous alteration environment and 
subsequent sample delivery. In the first model, the interior of a carbonaceous 
chondrite parent object is nearly homogeneously altered, and samples are derived 
following subsequent catastrophic break-up. In the second model, carbonaceous 
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chondrites are derived solely from regoliths, without sampling the deep interiors. 
For interior alteration, bounds on parent body size, composition, and heat source 
magnitude are given consistent with constraints on peak temperature and 
homogeneity. We present representative calculations showing the thermal effects 
of aqueous alteration reactions. Several mechanisms for introducing water into 
near-surface regions are illustrated for the regolith-alteration hypothesis, and we 
address the issue of long-term retention of ice on parent bodies. Finally, we 
discuss the remaining discrepancies between theory and observation, and we 
suggest revised models for the aqueous alteration environment of carbonaceous 
chondrites. 

ALTERATION IN CARBONACEOUS CHONDRITE PARENT BODIES 

Two important processes that affected early carbonaceous chondrite history 
were brecciation and aqueous alteration. Most, and possibly all, CM and Cl 
chondrites are breccias. As such textures have been interpreted as indicating that 
these meteorites have been exposed and reworked in a surficial regolith, it is then 
important to ask if aqueous alteration, too, occurred in a regolith environment 
[Kerridge and Bunch, 1979; Bunch and Chang, 1980]. Although many carbon- 
aceous chondrites are gas-rich [Goswami et ai, 1984], the "planetary” abundance 
pattern of noble gases does not support solar irradiation within a regolith. 

Petrographic observations of carbonaceous chondrites [ Richardson , 1978; 
Kerridge and Bunch, 1979; Tomeoka and Buseck, 1988] indicate that phyllosilicate 
formation and filling of fractures by mineralizing fluids occurred in several stages, 
consistent with the idea that mineralization and impact brecciation were 
contemporaneous. However, earlier alteration in asteroidal interiors cannot be 
ruled out, as some observations of compositional variations suggest that aqueous 
alteration of matrix materials may have largely preceded mechanical reworking. 
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McSween and Richardson [1977] found that matrix compositions in Cl chondrites 
are fairly constant within individual breccia clasts but vary between different clasts; 
similarly, host and clasts in CM chondrites also have distinct compositions 
[McSween, 1987]. Moreover, Richardson [1978] was able to describe the 
compositions of vein mineralizing solutions by local leaching of pre-existing matrix 
phyllosilicates, which suggests that filling of fractures postdated primary alteration. 

From the isotopic evolution of Sr in Cl carbonate, Macdougall et al. [1984] 
inferred that vein formation occurred very early, within 100 m.y. of the 4.55 b.y. 
formation time of unaltered chondrites. Fission-track compaction ages, 
representing the times at which altered phyllosilicate and isolated olivine grains 
came into final contact, range from 4.4 to 4.2 b.y. for CM chondrites [Macdougall 
and Kothari, 1976]; individual ages are accurate only to within 100 m.y. or so due to 
uncertainty in the initial 244p u /238u ratio. Taken together, these measurements 
provide an approximate upper limit of a few hundred million years for the duration 
of brecciation and alteration. 

Both brecciation and aqueous alteration of carbonaceous chondrites, then, 
must have occurred very early, either during or shortly after accretion, although the 
relative timing of these processes remains uncertain. Two broad classes of 
hypotheses that can account for aqueous alteration differ in the location of, and 
thus the heat sources and physical mechanisms involved in, the hydration of matrix 
minerals. In the first class, alteration may have occurred within planetesimal 
interiors, and brecciation produced during accretion into parent bodies, or 
alteration may have taken place within the interiors of the final parent bodies, with 
samples later excavated by catastrophic break-up. Scott and Rajan [1981] and 
Grimm [1985] have suggested similar models for ordinary chondrite metamorphic 
history. In a competing class of hypotheses, brecciation and alteration may have 
taken place contemporaneously within a regolith. Because these processes 
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occurred early in the history of the parent body, such a regolith could probably not 
be defined in the usual sense of a soil layer gradually excavated from underlying 
bedrock, but instead might include material above some depth to which the most 
recently accreted material was incompletely consolidated. These two classes of 
hypotheses may be termed interior alteration and regolith alteration, respectively. 

The extent of alteration within a particular environment warrants further 
consideration. Carbonaceous chondrite parent bodies may have experienced 
pervasive aqueous alteration under either interior alteration or regolith alteration, 
with the result that these zones were largely petrologically homogeneous, or 
alteration could have been limited, such that these regions contained significant 
amounts of unaltered materials. Two lines of evidence that bear on this problem 
are asteroid collision mechanics and carbonaceous chondrite chemistry and 
texture. 

Studies of asteroidal collisional histories [Daw's et a/., 1979; Housen et at., 
1979] suggest that catastrophic impacts, which result in complete fragmention of a 
body, may have occurred repeatedly for most asteroids over solar system history, 
although the largest objects (diameters from 300 to 1000 km) probably 
experienced few, if any, such events [Greenberg and Chapman, 1983]. Most 
sizable asteroids gravitationally reaccrete following catastrophic collisions, but 
some mass fraction is lost. Although surficial rocks may be accelerated to greater 
velocities and preferentially escape during catastrophic collisions, it is likely that 
some samples from the asteroidal interior would also be ejected, especially over 
several episodes of fragmentation and reassembly. Moreover, the object may 
eventually be disrupted and dispersed. From these considerations, we may 
reasonably infer that if a parent body has been totally shattered one or more times, 
then samples exist which show chemical and mineralogical features typical of 
conditions within the object's interior, if not representative of the entire range of 
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conditions. Strong heterogeneity in the extent of alteration in an interior-alteration 
scenario, then, seems unlikely on the grounds that it is unreasonable that only 
altered carbonaceous chondrites would be delivered from a shattered parent body 
without any samples of other kinds of interior mineralogy. A homogeneously 
altered parent body offers the most straightforward way to sample a parent body 
interior and obtain only altered materials. 

A regolith site for aqueous alteration may be excavated by either catastrophic 
or noncatastrophic collisions. Catastrophic events again raise the preceding 
question regarding sampling, leading to the identical conclusion that a 
homogeneous object is the most plausible condition. Alternatively, carbonaceous 
chondrites may be "chipped off" by noncatastrophic events that sample shallow 
depths of the parent object only (i.e., normal impact cratering). If the parent object 
experienced only these kinds of events over the age of the solar system, the 
implied strength against fragmentation probably is consistent only with a large 
asteroid. The regolith-alteration scenario, therefore, may be restricted to large 
objects, although the composition of the parent body's interior is unconstrained. 

The near one-to-one correspondence between petrologic type and chemical 
group among carbonaceous chondrites (i.e., the absence of unaltered and altered 
chondrites of the same group) argues further that processes which altered texture 
and mineralogy were pervasive, supporting either a homogeneous object or 
homogeneous regolith. In addition, although individual CM and Cl chondrites are 
largely brecciated, they almost invariably contain only altered materials. The only 
exceptions are a few C3 clasts in Murchison [Olsen et a/., 1988], a minor 
component that appears to have been a late addition from elsewhere after the bulk 
of Murchison was altered. 

The principal alternative models for the aqueous alteration environment under 
consideration here, then, are that of homogeneous alteration of an object's interior, 
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either in parent bodies or accreting planetesimals, or regolith alteration on a large 
object where the conditions at depth are unknown (Figure 6.1). Of course, the 
comparatively small number of samples, and the presumably smaller number of 
parent body impact events responsible for their excavation, make these 
hypotheses only the most likely out of a range of possibilities. 

CONSTRAINTS ON THERMAL MODELS 

Thermal models of ordinary chondrite parent bodies [e.g., Miyamoto et at., 1981] 
may be constrained to match peak temperatures inferred from the distribution of 
petrologic types among a particular chemical group, as well as cooling rates 
inferred from metallographic or isotopic closure methods. Among carbonaceous 
chondrites, differences in petrologic type may indicate variations in the extent of 
aqueous alteration rather than variations in temperature [McSween, 1979]. Peak 
temperature constraints are available for the carbonaceous chondrites, however, 
on a group basis. Moreover, hydrousranhydrous mineral ratios and isotopically- 
determined waterirock ratios provide additional constraints on prior history. 
Metallographic cooling rates, in contrast, are not a meaningful guide to parent body 
size or thermal history type for chondrites of petrologic type 3 or lower, because the 
low peak temperatures experienced by these meteorites lie within or below the 
metallographic closure interval [Miyamoto etal., 1981; Grimm, 1985]. 

Peak Temperatures 

The peak temperatures experienced by CM and Cl chondrites, while in contact 
with water, have been estimated by Clayton and Mayeda [1984]. The strong 
oxygen isotope fractionation between between calcite and hydrous matrix phases 
restricts the temperature during aqueous alteration to 25°C or less for CM 
chondrites. Under the assumption that Cl chondrites had an anhydrous precursor 
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similar to that observed in CMs, Clayton and Mayeda [1 984] derived a limiting 
alteration temperature of 150°C for Cl meteorites. These are the only available 
temperature estimates for the alteration process, although the higher temperatures 
inferred for Cl chondrites may be inconsistent with the observation that Cl 
phyllosilicates are more poorly crystallized than those in CM chondrites [Tomeoka 
and Buseck, 1988]. 

These bounds are applicable only for the period during which rock was in 
contact with water; if water was later expelled or exhausted by chemical reactions, 
subsequent temperatures would have been free to rise. Rietmeijer and Mackinnon 
[1985] inferred that temperatures of approximately 190°C for Cl and 260°C for CM 
chondrites were necessary to produce the observed ordering in poorly graphitized 
carbon. However, these estimates are accurate only to within 100°C, and 
temperatures may also be systematically shifted to lower values by catalytic activity 
or long heating times. An alternative peak temperature constraint is imposed by 
dehydration of matrix phyllosilicates. Textures characteristic of partial dehydration 
of serpentine-type phyllosilicates have been observed recently in two CM 
chondrites [Akai, 1988]. Although the search for similar textures in other 
carbonaceous chondrites is incomplete, it is clear that these meteorites did not 
experience temperatures significantly in excess of the dehydration interval, or else 
there would be no water-rich matrix. Akai [1988] experimentally reproduced the 
patchy, disordered and defect structures observed in two Antarctic CM chondrites 
by heating samples of Murchison to 500°C for one hour. Such textures were not 
observed at 400°C. The short duration of these experiments, however, makes 
extrapolation to geological time scales uncertain. Instead, we adopt a value of 
350°C as the value for the long-term dehydration temperature, based on measure- 
ments by Johannes [1968] of the equilibrium phase boundary for the reaction 
transforming serpentine and brucite to forsterite and water. We discuss this 
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reaction further in a later section. A temperature of 350°C may then be taken as the 
limiting value for Cl and CM chondrites following expulsion or exhaustion of H 2 0. 
For anhydrous CO chondrites, textural, mineralogical, and chemical studies limit 
the peak temperature to 450°C [McSween, 1 977]. 

Timing and Extent of Aqueous Alteration 

The strontium isotopic composition of vein carbonate in Cl chondrites 
[Macdougall et a!., 1984] and compaction ages for CM chondrites [Macdougall and 
Kothari, 1976] limit the alteration period to within a few hundred million years of 
accretion. Since we will argue below that the alteration reaction itself was rapid, 
this constraint applies principally to the duration of heat sources and the availability 
of water. It is assumed here that the entire matrices of Cl and CM chondrites are 
composed of altered materials; the required volume factions v h and v a of hydrous 
and anhydrous minerals are then in the ratios v h /v a ~ 1 .0 and ~ 0.5, respectively 
[Dodd, 1981]. CO chondrites are taken to be completely unaltered, so v h = 0. 

Water.Rock Ratio 

In addition to peak alteration temperatures, Clayton and Mayeda [1984] also 
estimated waterrrock oxygen mass ratios for aqueously altered carbonaceous 
chondrites. These ratios may be expressed as water volume fractions of 0.8 and ~ 
0.5 for Cl and CM chondrites, respectively. Such estimates correspond to the total 
volume of water flushed through a volume of rock over time, because no 
reasonable pore volume can hold 80% water. Therefore H 2 0 transport is indicated 
the Cl water:rock ratios. Because not all of the water that moves throughout the 
asteroidal interior will exchange oxygen isotopes with the rock, we may assume 
that these water volumes are lower limits. 
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FORMULATION OF THERMAL HISTORY MODELS 

Even within the restricted framework of the homogeneous interior alteration or 
regolith alteration models, uncertainty in the initial conditions generally prohibits a 
unique solution for all model parameters consistent with observational constraints. 
A senes of interior-alteration models is examined below to test for acceptable 
combinations of size, water volume fraction, and heat source abundance that 
satisfy peak temperature and homogeneity constraints. We consider only 
qualitatively the constraints of water:rock ratios. We complete the interior- 
alteration investigations with some simple tests of aqueous alteration chemical 
reactions on parent body thermal history. For regolith-alteration models, our 
approach is of necessity more exploratory. By explicit models, we seek to test the 
viability of possible mechanisms for introducing water into such an environment 
under representative conditions. 

In our basic model, carbonaceous chondrite parent bodies are taken to be 
initially cold mixtures of ice and anhydrous silicates heated by 26 AI and chondritic 
complements of long-lived radionuclides. When ice is melted, water may react with 
the rock to produce hydrous silicates. The effects of latent heat are included for 
both H 2 0 and silicates. Dehydration is assumed to occur at 350°C. Under 
sufficiently large temperature gradients, hydrothermal circulation may develop, 
redistributing heat so the interior reaches an isothermal state. Gas diffusion occurs 
as temperatures rise and increasing amounts of H 2 0 are vaporized. At high 
pressures, the overlying rock will fracture, and H 2 0 will be vented into the regolith 
or to space. Although we have attempted to identify several important processes 
that contribute to the thermal history of carbonaceous chondrite parent bodies, in 
practice not all elements of this model are evaluated simultaneously. 

While the mathematical details of the thermal history models are given in 
Appendix 6.A, supporting arguments for this model are presented in this section. 
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We discuss separately those processes occurring before the period of aqueous 
alteration and those taking place during alteration. Among the former processes 
are nebular condensation, accretion, and lithification of carbonaceous chondrite 
parent bodies. Among the issues important for the alteration stage are the heat 
sources responsible for volatile mobilization and mild metamorphism, the 
permeability of parent bodies undergoing hydrothermal circulation and gas 
transport, and the rate of aqueous alteration. 

Nebula Condensation 

The equilibrium condensation sequence [Grossman, 1972; Lewis, 1972; Goettel 
and Barshay, 1978] has provided a basis for understanding the radial 
compositional structure of the solar system. For potential aqueous alteration of 
carbonaceous chondrites, two crucial events were the formation of major 
H 2 0-bearing phases such as serpentine and talc at T < 400 K and the 
condensation of water ice at T < 200 K [Lewis, 1972]. The radial distances at which 
these phases appear may be estimated from models of the pressure and 
temperature structure of the solar nebula [Cameron and Pine, 1973; Wood and 
Morfill, 1988]. Lewis’s model shows the formation of talc in the region of the 
asteroid belt, while ice did not condense until nearly the distance of Saturn. By this 
model, hydrated silicates provided the original source of water in the region from 
Mars to Jupiter, including the putative asteroidal sources of meteorites. This 
assumption was implicit in the DuFresne and Anders [1962] model for 
carbonaceous chondrite parent bodies. 

Chemical and isotopic evidence in carbonaceous chondrites for the nebular 
equilibration of dust and gas has been critically reviewed by Bunch and Chang 
[1980], who suggested instead that direct condensation of ice and its accretion into 
parent bodies was possible. Recently, Prinn and Fegley [1988] have questioned 
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whether the equilibrium condensation assumption itself is valid, i.e., was chemical 
equilibrium reached within 10 5 -10 6 y, the lifetime of the solar nebula? By simple 
thermokinetic arguments, Prinn and Fegley showed that formation of hydrous 
minerals under nebular conditions would require time scales > 10 15 y for 
solid-solid reactions and 1 0 9 -1 0 10 y for gas-solid reactions. Even the more rapid 
case is 10 4 times longer than the lifetime of the solar nebula, providing a powerful 
argument that formation of primary hydrous minerals from the solar nebula at T < 
300 K is kinetically inhibited. On the basis of revised condensation and nebula 
models, Prinn and Fegley favored direct condensation of ice at the distance of the 
present-day outer asteroid belt. Considering the equilibration and thermokinetic 
arguments together, we assume that H 2 0 was incorporated into primary parent 
objects as ice instead of as an earlier generation of hydrous minerals. 

The initial proportion of ice and its degree of mixing with rock depends on the 
relative timing of ice condensation and parent body accretion. We consider a 
homogeneous accretion scenario, so that ice and rock are initimately mixed 
throughout the parent body. Dodd [1981] described carbonaceous chondrite 
matrix materials as "chemically unfractionated and the result of essentially 
complete condensation of nebular gas". It is not known, however, if this quenching 
extended to temperatures low enough for H 2 0 condensation, so that ice and rock 
would accrete together. Voids interpreted to have formed during evaporation of 
volatiles have been observed in the Orgeuil Cl chondrite [Tomeoka and Buseck, 
1988], but unambiguous textural evidence for the presence of ice-rock mixtures has 
not been found. An upper limit to the H 2 0 volume fraction of ~ 0.8 is given by 
cosmic abundances [Anders and Ebihara, 1982]; such high proportions of ice are 
appropriate for a cometary origin of carbonaceous chondrites at solar distances 
much further than the asteroids. At present, however, we limit the volume fraction 
of ice to 0.4 based on mechanical considerations (see below). 
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Accretion 

We assume that accretion is rapid compared with interior heating, so that the 
initial temperature profile is uniform and equal to the surface value in equilibrium 
with insolation (impact heating during accretion is negligible). Because plausible 
heat sources have durations of 10 6 - 10 7 y (see below), accretion times must be 
shorter than these values if this assumption is to be valid. Numerical simulations 
of planetesimal accretion at 2.7 A.U. [Greenberg et al., 1978] show that objects 
several hundred kilometers in diameter may form within ~ 10 5 y. Recent revisions 
to this algorithm [Patterson and Spaute, 1988] have resulted in a sharp upward 
revision of accretion times for the smallest objects, but the integrated time to 
accretion of large asteroids is relatively unaffected. On the other hand, analytical 
approaches [Safronov, 1979; Wetherill, 1980] predict ~ 10 7 y accumulation times 
for large asteroids. As an upper limit, asteroid accretion times are constrained by 
the growth time of Jupiter. By gravitational perturbations due to Jupiter itself, or by 
the scattering of Jupiter-zone planetesimals into the asteroid belt, encounter 
velocities between asteroids were increased such that further accretion was halted 
in favor of fragmentation and dispersal [Safronov, 1972, 1979]. Numerical studies 
by Davis et al. [1979] suggest that interactions with earth-sized Jupiter-scattered 
planetesimals raise typical asteroid encounter velocities to 5 km/s within 10 6 -10 7 y. 
These calculations are based on the presumed pre-existence of such 
planetesimals, which may themselves take ~ 10 8 y to accrete [Wetherill, 1980]; 
however, a runaway growth could occur here as well [Patterson and Spaute, 1988]. 
Since the time scale for accretion cannot be adequately resolved from the time 
scale for internal heating, the assumption of short accretion times must be viewed 
as ad hoc. Isotopic techniques cannot be considered to date accretion times 
reliably, in light of possible late disturbance by aqueous alteration. 

Following accretion, temperatures and compositional changes are assumed to 
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be the result only of the internal evolution of the parent object. There is some 
evidence among ordinary chondrites, however, to suggest that internal evolution 
was contemporaneous with further accretion [Scott and Rajan, 1981] or with 
collisional disruption and gravitational reassembly of the parent object [Grimm, 
1985]. We address contemporaneous alteration and accretion by finding the lower 
limit in size for objects that can satisfy the appropriate peak temperature 
constraints. 

Parent Body Compaction 

It is assumed that carbonaceous chondrites are fully lithified throughout the 
history represented by these calculations, although the mechanisms that could 
lead to consolidation and lithification in meteorites, such as overburden [Cain etai, 
1985], impact [Sugiura and Strangway, 1983; Sneyd et al., 1988], thermal 
recrystallization [Yomogida and Matsui, 1983, 1984], or pressure solution, are 
poorly understood. Most chondrites are coherent breccias in which clasts are also 
well lithified; the properties of these samples suggests that their parent bodies were 
well indurated. It is not known, however, if this induration was achieved before 
aqueous alteration. This assumption is critical, because if carbonaceous 
chondrites were still unconsolidated when H 2 0 was mobilized, differentiation of 
water and silicates may have occurred. 

An upper limit to the rate of upward expulsion of water from a rock skeleton may 
be specified if the silicates are completely deformable, so that pressures in the 
solid and fluid are the same. The separation velocity u is then given by Darcy's 
law, driven by a pressure proportional to the density difference Ap between the two 
phases: u ■ K Apg/rj, where K is the permeability, g is the gravitational 
acceleration, and rj is the viscosity of water. The Darcy segregation time to ~ R/u is 
the characteristic time to move water over a length scale R. Taking g to be the 
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average gravitational acceleration by volume within a uniform sphere of radius R 
and density p, g = it p G R, where G is the gravitational constant. With p = 3 Mg/m 3 
and t) = 10' 3 Pa-s, we find t D ~ 10' 13 / K, expressed in millions of years. If this 
quantity is small compared with the characteristic thermal diffusion time for a 
sphere, t ~ 0.1 R 2 /k, where k is the thermal diffusivity, then water-silicate separation 
may occur. Since k ~ 10’ 6 m 2 s- 1 , the thermal diffusion and Darcy segregation 
times are comparable for KR 2 ~ 10' 4 . At K = 10* 13 m 2 , then, upward expulsion of 
H 2 0 may begin at R ~ 30 km, with significant effects probably restricted to R > 100 
km. However, the permeability is poorly constrained (see below) and could vary 
by two orders of magnitude in either direction. A more realistic treatment of the 
problem of fluid expulsion from an elastic, deformable skeleton follows from soil 
consolidation theory [e.g., Lambe and Whitman, 1969] but at the expense of 
introducing additional unknown parameters. We find no reason to reject 
hypotheses that H 2 0 expulsion by silicate compaction is either fast or slow. Our 
assumption that it is slow is largely for convenience in the present model, although 
we have considered the effect of complete water-rock differentiation on 
carbonaceous chondrite thermal history in a preliminary report [ Grimm and 
McSween, 1988]; suggestions for future models with large water volumes are given 
in the Discussion below. The maximum water volume fraction in the fully-lithified 
case is then limited by the amount that will fit into the pores of a connected rock 
skeleton, approximately 40% by volume. 

Heat Source 

The heat source for primordial metamorphism and melting of meteorite parent 
bodies is still unresolved [Sonett and Reynolds , 1979]. However, two leading 
candidates have emerged: solar electromagnetic induction and the decay of the 
short-lived radionuclide 26 AI. 
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Electromagnetic heating of asteroids by the interaction of the early solar wind 
and solar magnetic field was first studied in detail by Sonett et al. [1970]. A 
scenario constructed by Herbert and Sonett [1979] for anhydrous objects shows 
marginal silicate melting on Vesta but none on Ceres; this is principally caused by 
the inverse square falloff of spherically symmetric plasma flow and associated 
induction heating. The thermal histories inferred from the observed 
spectrophotometric structure of the asteroid belt with solar distance [Gradie and 
Tedesco, 1982; Bell, 1986] may be a strong constraint on heating models. 
However, there are a number of unknown model parameters associated with the 
physical state of the early sun, and an asteroid's response is sensitive to its 
electrical conductivity. Therefore, the extent of electromagnetic heating and 
uniqueness of model scenarios are difficult to evaluate quantitatively. 

Aluminum-26 is the most favored of the potential short-lived radionuclides 
because of its expected high abundance, high decay energy, and reasonably long 
half-life [Sonett and Reynolds, 1979]. Since the discovery of its decay product 
(excess 26 Mg) in Allende by Lee et al. [1976], numerous thermal models for 
meteorite parent bodies have been constructed using 26 AI as a heat source (see 
Basaltic Volcanism Study Project [1 981] for a review). 

Because the 26 Mg excess was detected by Lee et al. [1 976] only in a refractory 
inclusion in Allende, the inferred initial abundance of 26 AI (initial 26 AI/ 27 AI = 
6x1 O' 5 ) may not be representative of the amount of 26 AI incorporated into parent 
bodies. In addition, it was immediately recognized by these workers that such 
levels of 26 AI would melt all but the smallest objects. Radiogenic 26 Mg had not 
been detected outside of refractory inclusions until Hutcheon et al. [1988] reported 
a 26 Mg excess in an anorthite-bearing chondrule in an ordinary chondrite, from 
which they inferred an initial 26 AI/2 7 AI of (8 ± 2) xIO' 6 . The discovery of radiogenic 
26 Mg outside of CAIs suggests that 26 AI was incorporated into parent bodies. 
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Although this new figure for initial 26 AI/ 27 AI may represent an upper limit, the value 
reported by Hutcheon et al. [1988] is not very different from that deduced from 
thermal models of ordinary chondrites. Under the assumption that the current 
collections of ordinary chondrites are representative samples of the H and L 
chondrite parent bodies, Miyamoto et al. [1981] found that an initial 26 AI/ 27 AI of 
5x10- 6 for both objects would have been sufficient to explain the observed 
distribution of petrologic types and their inferred peak temperatures. 

We believe that these modeling and experimental efforts support the case for 
26 AI as an important primordial heat source for chondrite parent bodies. The 
present model, then, is formulated in terms of radionuclide heat sources; however, 
the principal conclusions on the thermal evolution of carbonaceous chondrite 
parent bodies are not markedly affected by choice of heat source, because both 
26 AI and electromagnetic induction operate on time scales of ~ 10 6 y and the 
internal temperature profiles produced by both sources are broadly similar [Basaltic 
Volcanism Study Project, 1981]. 

Since the heating requirements for aqueous alteration of carbonaceous 
chondrites are modest, mechanisms that have been ruled out for the heating of 
meteorite parent bodies - in particular, impacts - could apply here. Lange et al. 
[1985] have experimentally demonstrated that shock dehydration of serpentine 
occurs over the pressure interval 25-59 GPa. Such release of H 2 0 could be 
responsible for later secondary alteration and veining, but is not relevent to the 
melting of primordial ice required for primary alteration. Shock melting of ice 
occurs over the interval 8-10 GPa; the corresponding minimum impact velocities, 
appropriate for a rock projectile and an ice target, are ~ 3 km/s [Chapman and 
McKinnon, 1986]. During the accretion period, encounter velocities are compar- 
able to the escape velocity of the largest object in the population [Safronov, 1972], 
which indicates that impact velocities between asteroids during accretion were only 
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several hundred meters per second at most. Impact melting of ice is therefore not 
likely during accretion. Encounter velocities among asteroids were gradually 
raised to the present value of ~ 5 km/s, however, by gravitational interactions with 
Jupiter or its scattered planetesimals. As described previously, consideration of the 
both the time scale for such planetesimal growth [Wetherill, 1980; Patterson and 
Spaute, 1988] and the time scale required to pump up asteroid velocities [Davis et 
al., 1979] only weakly constrains the time scale on which impact melting of ice can 
become effective to 10 6 - 10 8 y. Since accretion times are probably 10 5 - 10 7 y, it 
appears that impact melting can effectively melt ice in postaccretional regoliths. 
Carbonaceous chondrites show only weak evidence of shock, which suggests that 
they were not subjected to hypervelocity impact. Unfortunately, the correlation of 
shock facies [Dodd and Jarosewich, 1979; Sneyd et al., 1988] and shock pressure 
[Sears et al., 1984] established for ordinary chondrites cannot be used to estimate 
shock pressures for carbonaceous chondrites, because significant strain could be 
accommodated by ductile matrix or ice. The shock levels experienced by 
carbonaceous chondrites and the contribution of impact release of water to 
aqueous alteration is therefore an open question. Nevertheless, carbonaceous 
chondrites still require some additional heat source to explain mild static 
metamorphism of the CO group [McSween, 1977]. 

H 2 0 Transport 

When heating of the parent asteroid reaches several hundred kelvins, large 
internal redistributions of H 2 0 may occur due to evaporation of water in the deep 
interior, diffusion through the porous interior, and recondensation in a 
(comparatively) near-surface shell. In addition, under sufficient heating, pore water 
itself may become gravitationally unstable and undergo hydrothermal circulation. 
Both of these processes are sensitive to the permeability, a quantity which is not 
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reliably known. 

Sugiura et al. [1984] measured the permeability of Allende (CV3) to be * 10* 15 
m 2 , a rather low value characteristic of terrestrial silt or clay [Lambe and Whitman, 
1969] but consistent with the fine texture of the matrix. In the field, however, and 
presumably in asteroids, permeability is sensitive to the presence of fractures, 
which may greatly accelerate rates of fluid transport. Bunch and Chang [1 980] 
described CM chondrites as having textures similar to terrestrial hydrothermally 
altered tuffs and mafic rocks. A possible terrestrial analog to such rocks is the 
Waiora formation, a 600-800 m thick unit of pumice breccias and vitric tuffs in New 
Zealand's Wairakei hydrothermal region. Mercer et al. [1975] adopted K = 10' 13 
m 2 (at observed porosity 0.2) for the bulk permeability of this unit, which includes 
the contribution from fractures. The volcanic rocks of the Waiora formation have 
larger grain (pore) sizes than the Allende matrix, so the bulk permeability of the 
Waiora formation probably will exceed that of Allende's parent body, unless the 
size and density of fractures between these two brecciated units varies 
substantially. The value 10' 13 m 2 may be taken as a reasonable upper limit for the 
permeability of lithified asteroids. In the case of unlithified material, as discussed 
above, however, the permeability may have been as high as 10' 11 m 2 , an 
approximate upper limit for terrestrial sediments [Lambe and Whitman, 1969]. In 
addition, if the highly porous, silt-to-fine-sand texture of the lunar regolith [Costes 
and Mitchell, 1970] is any guide to asteroids, the near-surface permeability should 
lie in the range 10* 11 - 10* 12 m 2 . 

If fluid or gas pressures exceed the sum of the rock's tensile strength and the 
lithostatic pressure, fracturing will occur and H 2 0 will be vented. Compressive 
strengths of equilibrated ordinary chondrites measured by Buddhue [1942] range 
from 62 bar to 3.7 kbar. From a typical ratio of compressive to tensile strength of 10 
and typical tensile strengths of both sedimentary and crystalline rocks [Jaeger and 
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Cook, 1976], 100 bar (10 7 Pa) is taken as representative of the tensile strength of 
asteroidal interiors. Such a value implies that resistance to gas or fluid fracturing is 
dominated by the tensile strength (and therefore is approximately constant) for 
objects less than about 200 km in diameter; lithostatic pressure becomes 
increasingly important at larger sizes. 

Rate of Aqueous Alteration 

Aqueous alteration has profoundly affected the textures of carbonaceous 
chondrites, from millimeter-scale fractures to micron-scale matrix phyllosilicate 
grains. DuFresne and Anders [1962] suggested that water was available for > 10 3 
y, on the basis of the regularity of dolomite crystals in Cl chondrites. The actual 
aqueous alteration reactions in carbonaceous chondrites were complex, leading to 
hydrous minerals such as serpentines, septechlorites, and montmorillonites. In 
order to assess quantitatively the time scale for aqueous alteration, the reaction 

2 Mg 3 Si0 4 + 3 H 2 0 = Mg 3 Si 2 05(0H)4 + Mg(OH) 2 
Forsterite Serpentine Brucite 

is taken to be representative of carbonaceous chondrite aqueous alteration. Both 
hydration and dehydration rates for this reaction were studied by Wegner and Ernst 
[1983], who give the time for 99% hydration of olivine as a function of fluid 
pressure, temperature, and grain size. Even at the low fluid pressures relevant to 
this problem, hydration of mm-size grains at 0°C will be complete in ~ 10 4 y. For 
the micron-size grains of carbonaceous chondrite matrices, hydration should take 
only about a year at 25°C (an upper limit to the alteration temperature inferred from 
isotopic studies, discussed above). Olsen et al. [1988] also concluded that the 
duration of carbonaceous chondrite aqueous alteration was short, ~ 10 2 y. 
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Because aqueous alteration is fast compared with the thermal evolution time 
scales of asteroids, we may consider such reactions to be instantaneous in the 
context of thermal models. Although characterized by a different time constant, 
dehydration rates are also rapid [Wegner and Ernst, 1983], due to increased 
thermal activation. 

In these models we assume that the fluid phase is pure H 2 0. On CM and Cl 
parent bodies, mineralizing fluids had more complex chemistries, as evidenced by 
the fact that Ca-Mg-Fe carbonates and sulfates precipitated sequentially 
[Richardson, 1978]. Any kinetic effects of lowered activity of H 2 0 is ignored in 
these calculations. 

RESULTS FOR INTERIOR-ALTERATION MODELS 

Common parameters for both interior-alteration and regolith-alteration models 
are orbital semimajor axis 3 A.U., albedo 0.05, and emissivity 0.8 (estimated values 
for C-type asteroids). In order to obtain upper limits to hydrothermal convection, a 
reference permeability of 10* 13 m 2 at saturated 20% pore volume is used. 

Effect of H 2 0 on Peak Temperature and Limits on 26 Al Abundance 

The following series of calculations demonstrates the primary effects of water on 
thermal evolution. The pore volume is taken to be saturated with H 2 0; because 
this simplification easily leads to hydraulic fracturing at modest temperatures, 
venting of water is not included. Aqueous alteration and degassing are also 
neglected for purposes of these illustrations. The remaining variables are parent 
object diameter, water volume fraction, and initial 26 AI/2 7 AI ratio. In terms of these 
quantities, the problem may be formulated to find upper and lower limits to the 
radionuclide abundance at specified size and rock:ice ratio. The lower limit is that 
amount of 26 AI required just to melt water in the deep interior. Obviously, melting of 
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a small volume fraction of the object will not expose much of its rock mass to water 
or lead to a homogeneously altered interior, but such a level of 26 AI may be 
regarded as a lower bound. The upper limit to 26 AI is constrained by the peak 
temperature for each group described above. In order to estimate this bound 
conservatively for the homogeneous model, the volume-averaged peak 
temperature is used; core temperatures may be somewhat higher. 

The thermal history for an object 100 km in diameter and 30% ice by volume, 
which is brought just past the H 2 0 melting point by an initial 26 AI/ 27 AI ratio of 
1.6x10* 6 , is shown in Figure 6.2. Following a rapid (< 1 m.y.) rise in temperature to 
the melting point, the interior remains so for several million years while the latent 
heat of fusion of ice is absorbed. Under this minimum-melting scenario, only a 
negligable volume in the asteroid's center rises a few degrees above freezing and 
is altered. 

The thermal history of a similar object is given in Figure 6.3, but with an initial 
26 AI/ 27 AI ratio of 2.0x1 0* 6 . Following H 2 0 melting, the central temperature rises to 
310 K in about 5 m.y. Water is gravitationally stable and no circulation occurs. The 
volume-averaged peak temperature of all rocks exposed to liquid water is 298 K 
(25°C); these temperatures match the upper limit for CM chondrites. However, less 
than half of the volume has melted -- the outer 1 1 km stay frozen. This result 
illustrates a problem with constraining CM chondrites to small parent bodies: the 
bulk of the object must remain frozen in order to accommodate the narrow range of 
allowable temperatures, and the homogeneous alteration assumption is then 
violated. We return to this point below. 

Figure 6.4a shows a similar calculation applicable to Cl chondrites. An initial 
26 AI/ 27 AI ratio of 4.3x1 O' 6 results in melting of ice throughtout 70% of the object's 
volume, and the mean peak temperature is 423 K (150°C). Hydrothermal 
convection has been suppressed in this calculation. Figure 6.4b shows a 
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complementary model, also matching the Cl peak temperature (initial 26 AI/ 27 AI = 
6.0x10- 6 ), in which hydrothermal convection is allowed. Circulation of H 2 0 has a 
strong effect on the thermal evolution, homogenizing temperatures within the 
interior with the result that the object almost uniformly reaches the required peak 
temperature. Enhanced heat transfer allows melting of ice closer to the surface 
while suppressing strong temperature increases in the deep interior. Since the 
adopted approximation for convection yields an overestimate of the heat transfer 
(Appendix 6.A), Figures 6.4a and 6.4b bound the range of thermal histories in an 
homogeneous Cl parent body undergoing hydrothermal convection. 

The calculations above have been repeated for a range of rock:ice 
compositions and object sizes. Results are summarized in Figure 6.5. At a 
specified parent body diameter D, the lowermost curve gives the lower bound on 
26 AI concentration, i.e., the level required just to melt ice. The next higher curve 
gives the maximum initial 26 AI/ 27 AI which produces volume-averaged 
temperatures equal to the inferred CM peak temperatures, and the next 
corresponds to Cl. Hydrothermal circulation does not develop in CM models, but it 
does for Cl, and has been suppressed in this third curve as for the calculations in 
Figure 6.4a. The uppermost, dashed curve gives the upper bound to 26 AI 
abundance, given generous assumptions on the efficiency of hydrothermal 
convection in Cl objects. From Figure 6.5, it can be seen that more 26 AI is needed, 
of course, to attain the same peak temperatures at a higher initial ice volume 
fraction. The radionuclide abundance must also be larger at smaller diameters to 
compensate for the increased heat loss. At D = 1000 km, peak temperatures for all 
cases can be achieved by the decay of long-lived radionuclides alone, with the 
exception that a limited amount of 26 AI may be required in such a large Cl parent 
body if hydrothermal convection is particularly efficient. The interval between 
minimum- melting and CM peak temperatures is smallest at D = 100 km because 
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such objects are large enough to be insulating on the 26 AI decay time scale, but 
not on the decay time scale of K, Th, and U. 

A calculation for anhydrous CO chondrites has also been carried out for each 
parent body diameter in Figure 6.5. Since the observed porosity of carbonaceous 
chondrites is ~ 25% [Yomogida and Matsui, 1983], the initial void volume fraction is 
taken to be 0.25 to illustrate the difference in thermal evolution between dry and 
H 2 0-saturated cases. The similarity of maximum radionuclide abundances 
between Cl and CO chondrites suggests that water may play an important role in 
controlling their different thermal evolution; that is, water in Cl chondrites may have 
acted as a thermal buffer to prevent the metamorphism experienced by CO 
chondrites. Moreover, the required initial 26 AI/ 27 AI ratios (several parts per million) 
are within the observational upper limit established for ordinary chondrites by 
Hutcheon etal. [1988] (8x1 0* 6 ), and reasonable agreement with the value inferred 
by Miyamoto at al. [1981] (5x1 O' 6 ) from thermal modeling of ordinary chondrite 
parent bodies can be found near D=1 00 km. These arguments are more difficult to 
apply to CM chondrites, however, as their low inferred peak temperatures limit the 
initial 26 AI/ 27 AI ratios to a factor of 2-3 below Cl parent bodies for the same 
diameter. If CM parent objects are smaller, more 26 AI can be accommodated, but 
at the expense of a smaller altered fraction of the interior. 

These models explicitly satisfy peak temperature constraints only. It is possible 
also to consider, but only qualitatively, the waterrock ratio derived from oxygen 
isotope data. Because hydrothermal circulation does not develop in homogeneous 
CM objects, the bulk water content must be at least 50% by volume, a value 
beyond those considered here. This is a lower bound, because exchange of water 
between liquid and silicates during aqueous alteration has been neglected in 
these calculations. Because of the greater range in temperature for Cl chondrites, 
hydrothermal circulation is more likely, and hence it is possibile that the 
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time-integrated water:rock ratio for some samples could be large. 

Finally, the presumption of interior homogeneity must be tested. We assume 
that at least 50% of interior rocks must be exposed to liquid water, in order that 
delivered samples of altered carbonaceous chondrites are representative of the 
later-disrupted parent body. CM models satisfy this constraint at D = 300 km but do 
not at D = 100 km. For Cl models without hydrothermal circulation, D = 100 km may 
be accepted, but D = 30 km is marginal. At large diameters, melting extends closer 
to the surface for water-convecting interiors than nonconvecting ones. However, 
there is not much difference for smaller sizes, and so Cl models with hydrothermal 
circulation at D = 30 km still produce ~ 50% altered rocks. Therefore, in addition to 
the restrictions on initial 26 AI specified by Figure 6.5, homogeneous interior- 
alteration models may be constrained to parent bodies a few hundred kilometers in 
diameter and larger for CM chondrites and several tens of kilometers in diameter 
and larger for Cl objects. 

We described earlier a scenario in which interior alteration occurred in 
planetesimals, before incorporation into parent bodies. Such a model may be 
excluded for CM chondrites because of the large object diameters required for 
homogeneous alteration, although accretion of smaller Cl planetesimals plausible. 

Thermal Effect of Aqueous Alteration 

In the preceding models, we neglected H 2 0 chemical exchange in order to 
demonstrate the thermophysical effects of the presence of the liquid phase alone. 
Three calculations are given here in the context of interior-alteration models to 
illustrate the potential effects of chemical reactions. Since little gas diffusion occurs 
below 273 K, degassing is again neglected, but finite strength against gas fracture 
is included. Figure 6.6a shows temperature profiles in a minimum-melting scenario 
for a 100-km-diameter object heated at an initial 26 AI/ 27 AI ratio of 1.2x1 O’ 6 . The 
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object's porosity is 20% and is saturated with H 2 0. Melting of ice occurs over 
about a 3 m.y. interval. Since aqueous alteration is rapid, hydrous minerals are 
effectively formed instantaneously. The final hydrous:anhydrous mineral ratio is ~ 
1 (comparable to that observed in CM chondrites), and essentially all of the water 
at depth is exhausted by the reactions (Figure 6.6b). The temperature in the deep 
interior increases markedly (~ 150 K) due to the substantial heat of reaction 
released, but the peak temperatures are well below those at which dehydration or 
thermal recrystallization would occur. Because of the speed of the reaction, the 
peak temperature during alteration is taken to be ~ 0°C, although in reality heat 
release from progressive alteration might cause substantial temperature 
excursions while liquid water is still available. As a minimum-melting scenario, this 
CM model also suffers from a paucity of rocks exposed to water. However, 
because the oxygen isotope temperature constraints do not apply if water is rapidly 
removed by aqueous alteration, higher heat source concentrations would push 
outward the radius of the melting zone. 

As core temperatures reach 623 K (350°C), dehydration occurs. Figures 6.7a 
and 7b show temperature profiles at several times and the final composition profile 
at 10 m.y., respectively, for a similar object with an initial 26 AI/2 7 AI ratio of 5x1 O' 6 . 
Melting of ice, aqueous alteration, and initiation of dehydration all take place within 
1 m.y. The model parent body history becomes rather violent when dehydration 
reactions occur, because the pore volume is suddenly filled with near-critical H 2 0. 
Since the pressure in the deep interior is several kilobars, immediate H 2 0 venting 
of the interior ensues, guaranteeing that that rehydration cannot occur when 
temperatures again fall below 350°C. In the end, anhydrous minerals again 
dominate half of the volume, and the altered zone is restricted to depths between 1 
and 1 1 km. As discussed above, there is no evidence for widespread dehydration 
among CM chondites. Furthermore, oxygen isotope systematics preclude 
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rehydration at a temperature as high as 350°C, even if H 2 0 were not vented from 
the parent body interior. 

A similar calculation (D = 100 km, initial 26 AI/ 27 AI = 5x1 O' 6 ) for Cl chondrites is 
shown in Figure 6.8. The highest water volume fraction allowed (0.4) is needed to 
achieve near-complete alteration of the anhydrous rock. The heat of reaction 
released is sufficient to send the temperature to the dehydration point, but not 
through it. In this example, the dehydration reaction is a barrier to further 
temperature excursions. 

RESULTS FOR REGOLITH-ALTERATION MODELS 

The release of liquid water occurs readily in interior-alteration models because 
of the insulating effects of overlying rock. In a cold surface regolith, however, 
supplying liquid is more difficult. We have conducted several thermal history 
calculations for internally-heated parent bodies to illustrate the mechanisms by 
which liquid water may reach near-surface regions. The possible extent of impact 
melting within a regolith is briefly discussed. 

To what depth does regolith alteration extend? Housen et al. [1979] 
numerically modelled regolith formation on asteroids and concluded that typical 
regolith thicknesses are a few kilometers or less. However, such regoliths develop 
over the age of the solar system (or until the object is disrupted), a condition that 
does not apply to the early regolith in which carbonaceous chondrites may have 
been altered. An upper limit to the thickness of such an accretional regolith may be 
estimated by considering the depth of the crater resulting from the largest possible 
impact that does not disrupt the asteroid. This depth is derivable in principle from 
energy and momentum considerations [Greenberg et al., 1978; Chapman and 
McKinnon, 1986]. A simpler approach is just to use the largest crater depth 
observed for small bodies (e.g., Herschel on Mimas), which suggests that impact 
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turnover could extend to a depth of several kilometers. We therefore accept as 
equally likely any internal heating model which supplies water to within such a 
distance of the surface. Even if liquid water reaches only the base of the regolith, 
impact stirring may allow access to upper levels. 

Internal Heating 

Water may be introduced into the regolith by internal heating in three ways: by 
direct melting of local ice, which may lead to H 2 0 circulation and replenishment 
from below, by venting of liquid or vapor along fractures caused by failure under 
high pore pressure, or by vapor diffusion through existing pores and cracks. 
Following the arguments given earlier that regolith-derived carbonaceous 
chondrites require a large parent body, several models at D = 300 km are given 
which illustrate these mechanisms. Representative values in all cases for the initial 
26 AI/2 7 AI ratio, porosity, and water volume fraction are taken as 5x1 O' 6 , 0.25, and 
0.1, respectively. 

Figure 6.9a shows the thermal evolution of such an object in which 
hydrothermal circulation is suppressed. Melting of ice occurs deeper than about 3 
km. However, gas pressures lead to venting of H 2 0 after about 1 m.y., which 
rapidly results in 70% water loss from the object (Figure 6.9b). This water is 
removed from the calculation but is equivalent to a saturated surface layer 17 km 
thick. Therefore fractures initiating at depth can supply large quantities of steam to 
the regolith; condensation there would provide a low-temperature environment for 
aqueous alteration. Condensation must occur at least for the late veining, as Cl 
textures imply that volatiles responsible for fracture-filling must have been in the 
fluid phase, rather than as gases [S.M. Richardson, quoted in Kerridge and Bunch, 
1979]. 

Figure 6.10 illustrates a calculation identical to that of Figure 6.9, but with 
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hydrothermal convection. Beginning after 0.4 m.y., vigorous circulation and rapid 
heat losses hold down internal temperatures such that gas pressures achieve only 
minor venting. Melting of ice, however, extends even closer to the surface (1km) 
than in the previous example, well within our suggested accretional regolith. 

Lastly, an example of gas diffusion is given in Figure 6.11. This transport 
mechanism has been unimportant in models given heretofore because of 
sub-boiling liquid-saturated states, homogenization of H 2 0 distribution by 
hydrothermal convection, or exhaustion of water by chemical reactions. In the case 
shown, hydrothermal convection is suppressed so that the H 2 0 distribution is not 
homogenized. Under diffusion alone, water vapor is depleted from the deep 
interior and deposited in the outer 10 km, perhaps within the depth of impact 
gardening, in the short span of 1 m.y. Therefore, under different initial conditions, 
direct melting and circulation, venting, and vapor diffusion all appear to be capable 
of satisfying the requirements of regolith alteration powered by internal heating. 

Impact Heating 

As discussed above, impacts are insufficient to shock melt ice during accretion 
but become more effective when asteroid encounter velocities increase via 
gravitational interactions with Jupiter or its scattered planetesimals. Safronov 
[1978] and Kaula [1979] have studied impact heating of the Earth and Moon during 
accretion. A similar accretional-heating model for aqueous alteration of 
carbonaceous chondrites is conceptually simple, but the heat source has five free 
parameters (number of impactors as a function of diameter, velocity, and time; 
length scale for heat deposition and fraction of impact heat buried) as opposed to 
the single free parameter (initial 26 AI/ 27 AI ratio) required for radionuclide heating. 
Here we can at least test whether the impact-heating hypothesis may be rejected 
on the grounds that heat losses will be too rapid to allow aqueous alteration. The 
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characteristic cooling time fat a depth z is given by t ~ z 2 /k, where k is the thermal 
diffusivity. Given t £ 10 4 y for aqueous alteration (see above) and k ^ 10* 6 mZs* 1 (a 
typical value for solid rock) for regolith, then z < 1 km. At shallower depths, smaller 
values of for k are required to accomplish alteration, but there are no restrictions at 
greater depths. Using simple models for the distribution of impact heat with depth 
[Kaula, 1979; Bratt et ai, 1985], we find that that typical post-accretional impacts 
can easily melt ice to a depth of several kilometers. The actual cooling time will be 
affected by other contributions to heat flow, such as heat losses from regolith 
stirring and release of the heat of the alteration reaction, but these simple 
calculations indicate that post-accretional impact heating may allow aqueous 
alteration at depths of a few kilometers within carbonaceous chondrite parent 
objects. 

RETENTION OF H 2 0 IN CHONDRITE PARENT BODIES 

Since we have expanded upon suggestions that ice may have been the original 
form of H 2 0 incorporated into parent bodies [Bunch and Chang, 1980; Prinn and 
Fegley, 1 988], it is worthwhile to consider the long-term fate of H 2 0 in possible 
asteroidal sources. If aqueous alteration, degassing, or venting are not sufficient to 
deplete such objects of H 2 0 during primordial heating, then volatile loss occurs 
gradually over solar system history, and some fraction of the initial H 2 0 within 
meteorite parent bodies may still be present. 

Appendix 6.B presents a simplified treatment of the near-surface degassing 
problem, and gives a relation (equation 6.B7) for the present depth of H 2 0 loss as 
a function of permeability and orbital semimajor axis. It should be recognized that 
this permeability is the effective value, including any potential effects of free 
molecular flow, as discussed in Appendix 6. A. Solving equation (6.A17) for pore 
sizes of 1 mm and 1 pm, the transition temperatures to Knudsen flow are found to 
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be 250 K and 385 K, respectively. Since the latter value is significantly above the ~ 
160 K surface temperature at 3 A.U., free molecular flow can greatly enhance near- 
surface degassing. The depth of H 2 0 loss is plotted as a function of permeability 
for a pore size of 1 mm in Figure 6.12a, and for a pore size of 1 pm in Figure 6.12b. 
Because the pore size used is the smallest likely value, Figure 16.2b gives 
conservative upper limits to the depth of H 2 0 loss. 

At 3 A.U., degassing is extremely slow even under the most conservative 
assumptions, so any primordial ice should still be retained. These results contrast 
sharply with predictions by DuFresne and Anders [1962] that sublimation should 
limit the lifetime of ice on meteorite parent bodies to only a few million years. 
However, disturbances to the 87 Rb- 87 Sr system [Mittlefehldt and Wetherill, 1979] in 
some carbonaceous chondrites suggest that these meteorites were leached of Rb 
as late as 1 b.y. after their formation. The smooth 39 Ar- 4 °Ar release pattern in 
Murchison [Dominik and Jessberger, 1979] indicates that gas was gradually lost 
from this CM meteorite at low temperatures. Although primordial heat sources 
cannot operate on these time scales, these isotopic studies support the idea that 
volatiles were not rapidly depleted. 

Collisional evolution is an additional process that may result in more rapid H 2 0 
loss. Collisions will cause H 2 0 loss both by an increase in surface area due to 
comminution and by shock vaporization. The loss due to comminution is probably 
small, for two reasons. First, the observed size distributions in both asteroids 
[Veeder, 1986] and laboratory fragmentation experiments [Hartmann, 1969] 
indicate that mass is concentrated in larger objects, so that a typical volume is well 
shielded. Second, most large asteroids probably reaccrete within a few orbital 
periods of fragmentation, and perhaps within times as short as a few days [Grimm, 
1985]. Such short exposure times will further minimize the surface area available 
for H 2 0 loss. It is more difficult to evaluate the contribution of shock vaporization. 
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Using impedance-match calculations tabulated by Kieffer and Simonds [1980], 
typical present-day asteroidal collisions should result in stress levels of a few tens 
of GPa over a length scale comparable to the size of the smaller asteroid 
(projectile). Since impact vaporization of ice occurs at pressures less than 20 GPa 
[Chapman and McKinnon, 1986], it is apparent that the projectile can be completely 
devolatilized, as can two comparable-size bodies. 

These considerations imply that large, primitive asteroids, e.g., Ceres, are 
particularly attractive candidates for primordial ice retention because they have 
undergone fragmentation and reassembly at most a few times, if ever [Daws et ah, 
1979; Greenberg and Chapman, 1983]. Using Figure 6.12b as a guide, we 
suggest that such objects should still have ice presently accessible within tens of 
meters of the surface if their regoliths are poorly developed. Ice should 
nonetheless be retained below a well-stirred regolith, perhaps at a depth of 1 km 
[Housen et ai, 1979]. Slow loss of vapor from the interior might leave thin films of 
frost on the surface, as inferred from spectrophotometric observations of Ceres 
[Lebofsky et ai, 1981]. 

Of course, this hypothesis would best be tested by in situ sampling or a close 
flyby of candidate asteroids. Nevertheless, there are two potential experiments that 
may provide some constraints. Upper limits to the density of H 2 0 in the asteroid 
belt may be derivable from ultraviolet spectrometer measurements by the Voyager 
spacecraft and could provide a boundary condition on a combined model of H 2 0 
production and loss there. Furthermore, high- resolution imaging of Earth-crossing 
asteroids might directly resolve volatile loss and yield estimates of the H 2 0 
degassing rates of these objects near 1 A.U. 
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DISCUSSION 

The interior-alteration models given here illustrate the potentially important role 
of H 2 0 in the divergent thermal histories of Cl versus thermally metamorphosed CO 
and ordinary chondrites. Radiogenic heating can be buffered by the latent heat of 
fusion and large heat capacity of water, and by enhanced heat loss under 
hydrothermal convection. Therefore low-temperature aqueous alteration can occur 
instead of thermal recrystallization at the same heat source abundance. 
Hydrothermal convection may also be important in homogenizing interior 
temperatures and attaining locally high water:rock ratios inferred for Cl chondrites 
from oxygen isotope studies. 

The CM group is difficult to reconcile with the homogeneous interior-alteration 
hypothesis, however. At the same 26 AI levels inferred for other chondrites, the CM 
parent body must be small enough to lose heat rapidly in order to maintain low 
temperatures imposed by oxygen isotope data; in this case strong temperature 
gradients preclude a homogeneously altered interior. Because convection does 
not develop in the model CM objects, even under the generous assumptions used, 
there can be no local water:rock enhancement, and so the interior-alteration model 
might plausibly be rejected for this meteorite group. 

The present treatment of aqueous alteration is undoubtedly oversimplified. 
Since complete reaction occurs immediately, but only after total melting of ice, 
water is never free to circulate. Consequently, there is no opportunity to satisfy 
constraints on water:rock ratios. In reality, melting and alteration might be 
contemporaneous. Heterogeneity in alteration reactions and significant convection 
velocities could allow water to circulate without reacting (cf. terrestrial hydrothermal 
zones). 

If alteration occurred in regoliths, some additional inconsistencies may be 
qualitatively understood. Water can be introduced in great quantities from below 
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by venting or vapor diffusion, but probably at irregular intervals, and removal of 
reaction heat is more easily accomplished from a near-surface location. Alteration 
by local melting of ice still fails to satisfy observed watenrock ratios, unless 
hydrothermal circulation extends into the regolith. Some fluid-transport mechanism 
is required to supplement the original DuFresne and Anders [1962] model. 
Aqueous alteration by impact heating faces a similar problem, in that there is no 
H 2 0 enhancement over the bulk composition of the regolith and the projectile 
population, unless impactors become progressively more ice-rich with time, or else 
the parent body accreted sufficiently far from the sun that near-cosmic proportions 
of ice were condensed. 

All of the scenarios considered here assume that carbonaceous chondrite 
parent bodies accreted at asteroidal distances with ice proportions of a few tens of 
percent. Alternatively, the original parent objects could have resembled comets, 
and later evolved into asteroid-like bodies [Wasson and Wetherill, 1979]. Aqueous 
alteration and oxygen-isotope exchange could occur during water-silicate 
differentiation [Grimm and McSween, 1988], although an improved model might 
consider sinking of silicate grains through mobilized or partially melted ice. The 
mechanism for induration of cometary residues into lithified carbonaceous 
chondrites remains an open question, however. 

In spite of the relative detail of the preliminary models developed here, further 
improvements are necessary to satisfy the observational constraints provided by 
carbonaceous chondrites. Explicit calculation of oxygen-isotope exchange is 
needed to assess the total water volume present and the magnitude of variations in 
oxygen isotope abundance within the parent body. Better information on the 
rate-limiting processes for aqueous alteration is required - although alteration can 
occur rapidly, it may not always do so in nature. These efforts would lead to 
time-dependent modelling of aqueous alteration. On a more fundamental level, 
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criteria are needed to to distinguish the formation location (inner or outer solar 
system), alteration environment (regolith or interior), and lithification mechanism of 
carbonaceous chondrites. Accretion processes and heat sources are also 
incompletely resolved. 

CONCLUSIONS 

The following conclusions can be drawn from this study: 

(1) Mixtures of anhydrous silicates and ice are suitable starting compositions 
for carbonaceous chondrite parent bodies that have experienced aqueous 
alteration. 

(2) Aqueous alteration may have occurred within either asteroidal interiors or 
surficial regoliths. Chemical, textural, and mechanical arguments support a 
petrologically homogeneous interior for the interior-alteration model and a large 
parent body for the regolith-alteration model. 

(3) The fusion heat of ice, the high heat capacity of water, and the ability of 
circulating water to enhance rates of heat loss may all significantly contribute to 
thermal buffering of primordial heat sources for carbonaceous chondrite parent 
bodies. 

(4) The short duration of both plausible primordial heat sources and relevent 
chemical reactions easily allows aqueous alteration to occur within time constraints 
imposed by isotope systematics. In fact, exothermic hydration reactions probably 
require that water was gradually introduced into reaction sites to maintain low 
temperatures. Because the reverse reaction is endothermic, it is a thermal barrier, 
which may explain why little dehydration is observed in carbonaceous chondrites. 

(5) For the interior-alteration model, bounds on initial 26 AI/ 27 AI ratios for Cl 
chondrites are found comparable to those inferred or observed for both CO and 
ordinary chondrites. Model Cl objects must be greater than several tens of 
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kilometers in diameter to satisfy the interior homogeneity condition. CM chondrites 
require lower 26 AI/ 27 AI ratios, and model objects of this group do not develop the 
hydrothermal circulation that is probably necessary to satisfy water:rock ratios 
implied by oxygen isotope data. 

(6) Regolith alteration may be driven by impacts or by internally-generated heat 
sources. In internally-heated models, temperatures at depth may reach higher 
levels than allowed for interior-alteration models. Large quantities of water may be 
supplied to the regolith from below by hydrothermal circulation, vapor diffusion, or 
venting. Impacts can provide sufficient heat after parent bodies have accreted but 
may not be able to deliver the required water volumes. 

(7) Retention of primordial ice in carbonaceous chondrite parent bodies is 
probably not limited by sublimation rates or by collisional comminution, but by 
shock vaporization. 

(8) The present models cannot quantitatively account for all observational 
constraints imposed by different carbonaceous chondrite groups. 
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Appendix 6.A: 

Mathematical Formulation of the Carbonaceous Chondrite Thermal Model 


The mathematical and numerical details are given here of the model used to 
calculate the internal evolution of carbonaceous chondrite parent objects. The 
principal assumptions adopted to arrive at this model are described in the text; only 
supporting details are given below. A summary of nomenclature is given at the end 
of this appendix, and all units are MKS unless noted otherwise. 


Heat Flow 

The spherically symmetric heat-conduction equation with internal heating 


|I i 3 (r 2 K fE. )+ J3 

3t r 2 3r 3r p C 


(6.A1) 


is solved by an implicit finite-difference method. The quantities in equation (6.A1) 
are temperature T, radius r, thermal diffusivity k, density p, specific heat c p , and 
heat production per unit mass Q. The thermal conductivity is given by k = p c p k. 
The parameters that control thermal diffusion (k, p, and c p ) are taken as simple 
weighted sums of the contributions of rock, gas, and ice or water. The natural 
weights for k and p are relative volume, whereas the relative mass contributions 
govern c p : 


k = 



X 


v j 



c p = 


X c p m j 

J V J J 


X 

J 


m. 


(6.A2) 
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where Vj and rrij are the volume and mass fractions, respectively, of the j th 
component. From equations (6.A2) follows 


K 5= 



£ p jVi 


(6.A3) 


for the effective thermal diffusivity. More precise formulations of effective thermal 
conductivity in two-phase media exist, such as the Maxwell formula [Carslaw and 
Jaeger, 1947], but the differences between (6.A2) and these formulations are minor 
for the present purpose. Temperature-dependent formulae for the thermodynamic 
and transport properties of rock and H 2 0 are given in a later section in this 
appendix. 

The heat production Q is the sum of the heat production of individual radioactive 
elements, weighted by their time-dependent mean abundances in carbonaceous 
chondrites: 


Q= mr lA oV V (6.A4) 

j i i 

where A 0 j is the initial abundance, Q 0 j is the initial heat production per unit mass, 
and Xj is the decay constant of the j th radionuclide, and t is time. Since H 2 0 is also 
present, an overall weighting factor m r , the mass fraction of rock, is required. The 
abundance of 26 AI, the primary heat-producing element in the model, is computed 
from the product of the assumed initial 26 AI/ 27 AI ratio and the observed mean Al 
content of CM chondrites. The relevant elemental abundances are given in Table 
6.A1. 
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The surface temperature is specified by solving the energy-balance equation 


-^(l-A) = eoT 4 +k^l I=R (6.A5) 

4d 


where S 0 is the solar constant, d is the orbital semimajor axis, A is the albedo, e is 
the emissivity, a is the Stefan-Boltzmann constant, and R is the object's radius. The 
loss of vapor is assumed to occur slowly, and its latent heat is neglected in (6.A5). 
This assumption follows from consideration of the characteristic time scale for 
surface degassing (Appendix 6.B) and has been verified by numerical exper- 
iments. 


Ice - Water Phase Transition 

Melting of ice is assumed to occur at a temperature of 273 K, independent of 
pressure. Since lithostatic pressures exceed a kilobar only for the largest 
asteroids, the error in this approximation is generally less than 10 K. During the 
phase transition, thermal energy is absorbed as the latent heat of fusion: 


AT = 


MjLf 



(6.A6) 


where AT is the temperature rise suppressed during melting, L t is the latent heat of 
fusion of ice, and M and Mj are the total mass and mass of ice, respectively. Upon 
freezing, AT is the suppressed decrease in temperature. Temperature changes in 
regions of the phase transition are tracked using the algorithm of Reynolds et al. 
[1966], which uses a partial melt state and associated partial release and absor- 
ption of fusion heat to avoid numerical instabilities associated with finite time steps. 
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Hydrothermal Convection 

Closed circulation of heated water may occur within the asteroid under favor- 
able conditions and provide significant heat transfer. A simple parameterized 
convection scheme is adopted, similar to calculations of the thermal histories of the 
Earth [Sharpe and Peltier, 1978], the Moon and the terrestrial planets [Cassen et 
ai, 1979; Schubert et al., 1979; Turcotte et ai, 1979], and the Galilean satellites 
[Thurber et ai, 1980]. 

Buoyant instability in a porous layer of thickness r w heated from below, with 
permeability K and average thermal conductivity k, occurs when the Rayleigh 
number 


KATr 


Ra 


w ° w p 


w 


w 


(6.A7) 




exceeds the critical value Ra cr = 4k 2 = 39 [e.g., Turcotte and Schubert, 1982], 
where p w is the reference density of water, g is the gravitational acceleration, a w is 
the thermal expansion coefficient of water, AT is the temperature drop across the 
layer, and ri w is the viscosity of water. The efficiency of convection is measured by 
the Nusselt number Nu, which is ratio of the total heat transported out of the 
convecting layer to the heat transported by conduction alone, Nu = q/q c . The 
relation between Ra and Nu has been numerically investigated for 
two-dimensional hydrothermal convection by Straus [ 1974]. These results are 
approximately described by the emprical formulae 


Nu=Ra' 


1.3 


Nu = 1.6 Ra' 


0.6 


Ra' <2 
Ra' > 2 


(6.A8) 
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where Ra' = Ra / When the layer is heated from within instead of from below, 
the reciprocal of the nondimensional temperature difference across the layer 
measures steady-state convective efficiency [Turcotte et al., 1979]; however, Nu 
may be used for the transient problem when the heat production per time step is 
small [Cassen et al., 1979]. 

In the parameterized convection approach, an equivalent conductive solution is 
sought that has the same heat transfer as the convective system; this may be 
accomplished by solving (6.A1) with an effective thermal diffusivity Nu k [Sharpe 
and Peltier, 1978]. The ’actual' temperatures are not those of the equivalent 
conductive solution, however, but must follow an adiabat within the convecting 
region. Since the slope of the H 2 0 adiabat is small at the low gravitational 
accelerations relevant to this problem, the circulating region is approximately 
isothermal; its temperature is calculated by matching the total thermal energy from 
the equivalent conductive solution [Cassen et al., 1979]. A boundary layer of 
thickness ~r w /A/t/ exists at the top of the convecting layer and is neglected in 
computing ’actual’ temperature profiles; this approximation is valid for Nu » 1. 

Since convection is not limited to upper layers in the object, the effects of 
spherical geometry will strongly influence the growth of density perturbations and 
the resulting flow. An analysis of the stability of an internally heated porous sphere 
has not been found in the literature. However, comparison of the conditions of the 
onset of convection [Chandrasekhar, 1961] for an internally heated fluid sphere 
and for a fluid layer heated from below (with mean gravitational acceleration of the 
sphere) shows that, for the same temperature drop, convection in the sphere 
requires Ra' approximately twice as large as for the layer. The initiation of 
instability in the sphere is inhibited by the smaller buoyancy forces resulting from 
decreased relative volume and lower gravity in the deep interior. Furthermore, 
Turcotte et al. [1 979] showed that steady-state heat transfer is more efficient for the 
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layer than for the sphere at the same value of Ra'. Use of (6.A7) - (6.A8), therefore, 
must also yield an overestimate of the convective efficiency of the porous sphere. 
The error is conservative, however, because an upper limit on the efficiency of 
heat transfer must be used to constrain upper bounds on the size and heat 
production of carbonaceous chondrite parent objects. The layer approximation 
sati$fies this constraint and renders the problem easily tractable. 

A further limitation to the parameterized convection approach is that the stability 
analysis which leads to the derivation of Ra requires that all properties are 
constant, with the exception of a linear temperature-dependent variation of density 
(Boussinesq approximation). Such an approximation is not always satisfied in 
hydrothermal systems with temperature variations of several hundred degrees; 
density variations are of order unity, and viscosity variations may span an order of 
magnitude. Terrestrial hydrothermal convection in a medium of variable 
thermodynamic properties was studied by Straus and Schubert [1977], who found 
that Ra cr may be much smaller than the Boussinesq value (~39). These results 
may be applied to asteroidal interiors by noting from (6.A7) that the thickness of a 
layer for which instability occurs scales inversely with gravity. Employing the 
scaled layer thickness and full temperature interval over which H 2 0 is liquid, we 
find from the work of Straus and Schubert that Ra cr may be significantly reduced for 
objects with diameters greater than 300 km. By using volume-averaged properties 
in (6.A7), some of these effects may be suppressed; nevertheless, this nonlinearity 
remains a potential source of error since it does not lead to upper limits to heat 
production. 

Hydrothermal convection will be limited by the degree of H 2 0 saturation of the 
pore volume. If water does not completely fill the pores, the medium may be treated 
as one of reduced effective permeability. The effective porosity in the unsaturated 
state is just the water volume fraction v w , and, using the Kozeny-Carman relation 
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[e.g., Lambe and Whitman, 1969], the effective permeability may be expressed as 


K=K 0 


e/(l+e) 

e 0 /(He 0 ) 


(6.A9) 


where e is the void ratio ( e = <j) / (1 -<>), where <> is the porosity) and the zero 
subscript denotes the saturated state. Since the variation in Kq between different 
rock types spans several orders of magnitude, the onset of convection (and the rate 
of H 2 0 vapor loss, discussed below) is very sensitive to the permeability. 
Arguments for the choices of this parameter are discussed in the text. 

The key assumption in evaluating effective permeability in unsaturated states is 
that water remains continuously interconnected within the pores. In soil 
mechanics, the maximum capillary head defines the height over which water 
remains in contact following draining of a soil [Lambe and Whitman, 1969]. Since 
capillary head is inversely proportional to both gravity and pore dimension 
[Batchelor, 1967], very large capillary zones, in the range 1 - 1000 km for pore 
dimensions of 1 pm to 1 mm, may exist in carbonaceous chondrite parent objects 
(g ~ 10‘ 2 m s' 2 ). Even if water were not continuously interconnected, it is 
conservative to assume so, for reasons given previously, in order to provide 
favorable conditions for convective heat transfer. 


Gas Diffusion 

The principal assumption of the vapor transport model is that the interior is in a 
state of quasi-static equilibrium. Under these conditions, the gas pressure always 
adjusts rapidly to the local temperature and H 2 0 volume, which change in 
response to diffusion of gas and thermal energy. 

From the Clausius-Clapeyron relation, the vapor pressure has the approximate 
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P =Pe (6. A 10) 

vap o v 7 

where P 0 and T 0 have dimensions of pressure and temperature, respectively. By 
fitting data given in Keyes [1928] and Washburn [1928], we adopt P 0 = 3.58x1 0 1 2 
Pa, T 0 = 6140 K for the vapor pressure over ice, and P 0 = 4.70x1 0 10 Pa, T 0 = 4960 
K for the vapor pressure over water. 

The equation of state used for H 2 0 vapor or gas is 

PV = nRT ( T < 500 K) 

2 < AU > 

(P + ^)(V - nb) = nRT ( T > 500 K) 

V 

where P is the ambient vapor or gas pressure, V is volume, n is the number of 
moles of H 2 0, R is the universal gas constant, and a and b are the van der Walls 
coefficients. From (6.A1 1) the gas density p g = n Ml V may be found (M is the molar 
mass); vapor-condensate mass balance then yields the mass of gas M g in a 
specified volume of the asteroidal interior; 


M = pV 
g K g g 


P C V<D-M 0 

Pc 

— 1 


(6.A12) 


where V g is the volume of gas, p c is the density of condensate (ice or water), V is 
the total volume, <|> is the porosity, and M 0 is the total mass of H 2 0. 
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The diffusion of H 2 0 through the interior must satisfy the equation of continuity 


T^(r 2 pu)=-fr (6 - A13) 


and Darcy’s law 


K 9P 
H 9r 


(6.A14) 


where u is the vertical velocity, rj is the viscosity, K is the permeability, and 
gravitational forces, which are small, have been neglected. Equations (6.A13) and 
(6.A14) may be combined into the nonlinear equation 


3p = J_^_/ 2 Kp_dP\ 

at r 2 0r V n 3r' 


(6.A15) 


which is solved by an explicit finite-difference method. The gas mass transferred in 
a given time step is AM g = Ap g V g . An advective temperature change 

AM L v 

AT = -— S-L (6.A16) 

Mc pg 

is caused by the net evaporation or condensation of vapor at a particular location, 
where L v is the latent heat of vaporization. Equation (6.A16) holds only when 
condensate is present and the small change associated with the product of the 
gas’s heat capacity c pg and the temperature difference between adjacent grid 
locations is neglected. 
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In this treatment of the gas diffusion problem, H 2 0 gas and condensate are 
considered immiscible fluids. The permeability is calculated from (6.A9), using v g 
for the effective porosity and 1-v r for the reference porosity. At low temperatures, 
however, the efficiency of diffusion may be greatly enhanced by free molecular 
(Knudsen) flow, which occurs when the mean free path of a gas molecule is 
significantly larger than the characteristic pore dimension: 

— »d (6.A17) 

where is Avogadro's number, b is the radius of the H 2 0 molecule, and d is the 
pore size. Equation (6.A17) may be solved as an equality for the approximate 
transition temperature to Knudsen flow under equilibrium conditions. Following 
Sugiura et al. [1984], free molecular flow may be parameterized as an increase in 
the effective permeability inversely proportional to the vapor pressure at 
temperatures below the transition temperature. In practice, Knudsen flow is of little 
importance during the active thermal history of the asteroid but becomes relevant to 
the treatment of secular degassing (Appendix 6.B). 

Venting of H 2 0 

Gas diffusion obeying (6.A14) will generally be smooth, unless pressures 
become so large that failure occurs. In this case H 2 0 may move rapidly to upper 
levels in the asteroid where it may recondense in available pores or escape 
directly into space. Failure occurs when the pore pressure exceeds the sum of the 
lithostatic pressure and the tensile strength of the rock [Jaeger and Cook, 1976]: 


P > P| ith + x 


(6.A18) 
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The lithostatic pressure is given by 


R 

P lith = f ,tG Jpr rdr (6 - A19) 

When gas fracturing occurs, H 2 0 has a direct path to the surface and might be 
free to boil away completely in an isothermal environment. However, advective 
heat losses will limit the rate of venting. In practice, (6.A16) is solved for the vented 
gas mass due to a temperature increase in a given time step, and the temperature 
is then reset to its original value, so that heating at a given location is balanced by 
venting of gas. As long as condensate is available, the temperature is constant 
and the pressure is maintained at the failure strength. The vented H 2 0 may be 
expressed as the thickness of a saturated surface layer of equal volume as a 
convenient measure of the depth to which significant liquid water may be available 
in the regolith. 

An additional failure mode, hydraulic fracturing, is possible when thermal 
expansion would lead to a volume of water in excess of the available pore volume. 
In this case, the pressure is given by 



(6.A20) 


where B is the bulk modulus of water, <|>V is the total pore volume, and V w is the 
volume to which water would expand. If (6.A1 8) is satisfied, the difference between 
these volumes is lost by hydraulic fracturing. 
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Hydration and Dehydration Reactions 

The formation of serpentine and brucite from forsterite and water is taken to rep- 
resent the reactions characterizing aqueous alteration of carbonaceous chondrites: 

2 Mg 2 Si0 4 + 3 H 2 0 = Mg 3 Si 2 0 5 (0H) 2 + Mg(OH) 2 (6.A21) 

The water:rock volume ratio of the reactants, which is denoted p, is 0.63, and the 
volume ratio of the products to the reacting forsterite, which is denoted y, is 1.56. 
Assuming that net volume changes between the silicate phases are accom- 
modated entirely by the pore volume, i.e., there is no net expansion of the asteroid, 
the volume balance for (6.A21) is: 

Av a + Avj, + Av w + Av v = 0 Av w /Av a = (3 Avh/Av a = -y (6.A22) 

where Av is the change in volume fraction and the subscripts a, h, w, and v imply 
anyhdrous minerals (forsterite), hydrous minerals (serpentine and brucite), water, 
and voids, respectively. Note that <|> = v w + v v . For the hydration (retrograde) 
reaction, two cases may be distinguished. If v w /p < v a , then incomplete alteration of 
anhydrous minerals will occur, and all water will be exhausted. If v w /p > v a , then 
complete alteration of all anhydrous minerals occurs, and some excess water 
remains. Using only the subscripts, the changes in volume fraction for incomplete 
hydration are 


Aw = -w Aa = -w/p Ah = yw/p Av = -w(l+p-y)/p (6.A23) 


and the fractional volume changes for complete hydration are 
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Aa = -a Aw = -(3a Ah = 7 a Av = -a(y-p-l) (6.A24) 

From the standard heats of formation tabulated by Robie and Waldbaum [1968], 
the heat of reaction released at 25°C is H = 233 kj/kg. The increase in temperature 
is 

AM, H 

AT = — — — (6.A25) 

Me 

P 

where AMh is the mass of hydrous minerals formed by the reaction. Since the 
reaction time is small compared with the thermal diffusion time (see text), formation 
of hydrous minerals and release of reaction heat are assumed to occur 
instantaneously upon complete melting of ice. 

The equilibrium phase boundary for this reaction has ben determined by 
Johannes [1968] to lie at ~360°C at a fluid pressure of 1 kbar. Fluid pressures in 
asteroids will usually lie below this value, even if they approach lithostatic 
pressure. Therefore, 350°C is taken to be the temperature of the dehydration 
(prograde) reaction. All hydrous minerals are assumed to dehydrate 
instantaneously, resulting in the fractional volume changes 

Ah = -h Aa = h/y Aw = (3h/y Av = h(y-p-l)/Y (6.A26) 

The changes in the enthalpies of forsterite and brucite at 600 K are negligible, and 
the enthalpy of water changes by only 14%. The high-temperature enthalpy of 
serpentine was not given by Robie and Waldbaum [1968], but changes in other 
phyllosilicate minerals are also negligible. The heat of the dehydration reaction at 
350°C, then, is assumed to be -233 kJ/kg. When the temperature falls below 
350°C, rehydration of the silicates occurs. In practice, temperature changes are 
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evaluated by the same algorithm used for melting and freezing of ice, by using a 
partial reaction state to avoid numerical instabilities. 

Temperature-dependent material properties 

The low-temperature thermal diffusivity of rock is taken to be that measured by 
Yomogida and Matsui [1 983] for Leoville (CV3): 

^-3.02x10-7 + 2.78x10-4/7 

Since the mineralogies of unaltered and altered rock are considered to be roughly 
those of olivine and serpentine, the density of rock p r is taken to be either 3.3 
Mg/m 3 or 2.5 Mg/m 3 , respectively. Specific heat measurements of carbonaceous 
chondrites are lacking, so c pr is taken to be 700 J/kg-K, a typical value for ordinary 
chondrites [Alexeyeva, 1 958; Yomogida and Matsui, 1 983]. 

Since 400 K was the highest temperature at which measurements on Leoville 
were performed by Yomogida and Matsui [ 1983], k r is dominated by lattice phonon 
conductivity. At higher temperatures, radiation becomes important. On the basis of 
thermal measurements of olivine and pyroxene, Schatz and Simmons [1972] 
include a radiative conductivity term that is directly proportional to temperature 
above 500 K. Assuming that the total thermal conductivity between 500 K and 
1000 K of the minerals measured by Schatz and Simmons [1972] and Leoville 
differ by a multiplicative constant, the total (lattice + radiative) conductivity k r ' is 

k r ’ = k r + 2.1x1 O' 3 (T- 500 K) 
at temperatures greater than 500 K. 

The thermal conductivity and heat capacity of hexagonal ice are given by 
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Hobbs [1974]: 

kj = 0.465 + 488 /T 
c pi = 7.67 T 

and, for simplicity, the density of ice is taken to be pj = 1.0 Mg/m 3 . Since 
amorphous ice is unstable at T > 135 K [Klinger, 1980], its different thermal 
properties need not be considered at temperatures typical of asteroidal surfaces. 

A fourth-order polynomial was constructed for the density of water, by fitting 
measurements summarized by Keyes [1928]: 

p w = -221 + 13.1 T - 5.07x1 O’ 2 T 2 + 8.49x1 0* T 3 -5.48x1 O' 3 T 4 

Gaseous H 2 0 is taken to behave ideally below 500 K and as a van der Waals gas 
at higher temperatures; the density p g may be calculated from the relevant equation 
of state. Details are given above in the section on gas diffusion. 

The thermal conductivity of H 2 0 liquid and gas are given by Touloukian et al. 
[1970]: 

kw = -0.581 + 6.34x1 0*3 T - 7.93x1 O' 3 T2 T < 41 0 K 

kw = -0.1 42 + 4.1 2x1 O' 3 T - 5.01 xl O' 3 T2 T > 41 0 K 

kg = -1.43x10*2+ 1.02x10- 4 T 

The specific heat of water is given by Touloukian and Makita [1970], and the 
specific heat of gaseous H 2 0 was taken from Zemansky and Dittman [1 981]: 
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C pw = 4200 T < 41 0 K 

c pw = -4.67x1 04 + 333 T - 0.731 T2 + 5.40x1 0*4 T3 T > 41 0 K 

Cpg= 1680 + 0.552 T 

where the small variation in the specific heat of water below 400 K and 
higher-order terms in the specific heat of H 2 0 gas have been neglected. 

The viscosity of H 2 0 liquid and gas is given by Touloukian et al. [1975]: 


T| w = 2.414x10-5 1 0247.8/(T-140) 
Tig = 8.04x1 0-6 + 4.07x1 0-8 T 


Nomenclature 

a first Van der Waals coefficient, 0.546 m 6 mol' 2 
a (subscript) anhydrous rock 

A radioactive element abundance, wt. fraction 

A albedo 

b second Van der Walls coefficient, 3.05x1 0 -5 m 3 mol* 1 
b radius of water molecule, ~ 5x1 0- 1 1 m 

B bulk modulus of water, 2x1 0 9 Pa 

c (subscript) H 2 0 condensate (water or ice) 

c p specific capacity at constant pressure, J kg* 1 K' 1 

d orbital semimajor axis, m 

d grain or pore diameter, m 

D object diameter, m 

e void ratio 

g gravitational acceleration, m s’ 2 

g (subscript) gas 

G gravitational constant, 6.67x1 O' 11 Nm 2 kg- 2 
h (subscript) hydrated rock 
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H heat of reaction for 2 Fo + 3 H 2 0 = serp + bru, 2.33x1 0 5 J kg* 1 
i (subscript) ice 

j summation index for bulk rock-H 2 0 composition or radionuclides 

k thermal conductivity, W/m-K 

K permeability, m 2 

Kq reference permeability, m 2 

L f latent heat of fusion of H 2 0, 3.34x1 0 5 J kg* 1 

L v latent heat of vaporization of H 2 0, 2.26x1 0 6 J kg* 1 

m mass fraction 

M mass, kg 

M molar mass of H 2 0, 1 8x1 0* 3 kg mol* 1 

n moles 

N a Avogadro’s number, 6.02x1 0 23 molecules mol* 1 
Nu Nusselt number 

P pressure, Pa 

P 0 pressure parameter for vapor pressure equation 
Q radionuclide heat generation, W kg* 1 

r radius, m 

r (subscript) rock 

r w radius of liquid water, i.e., width of potential convecting layer 
R object radius, m 

R universal gas constant, 8.314 J mol* 1 K* 1 

Ra Rayleigh number 

Rei Rayleigh number, expressed as multiple of critical value 
S 0 solar constant, 1 360 W m* 2 

t time, s 

T temperature, K 

T 0 temperature parameter for vapor pressure equation 
u velocity, m/s 

v volume fraction 

v (subscript) void 

V volume, m 3 

a volumetric coefficient of thermal expansion, K* 1 
p water/rock volume fraction for hydration 

y volume expansion factor of hydrous minerals 

e emissivity 
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rj viscosity, Pa-s 

k thermal diffusivity, m 2 s _1 

X radionuclide decay constant 

p density, kg/m 3 

a Stefan-Boltzmann constant, 5.67x1 0 8 W nr 2 K* 4 

x tensile strength, Pa 

<|> porosity 
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Appendix 6.B: 

Degassing and Secular Retention of H 2 0 


The model derived in Appendix 6.A may be used to calculate the diffusion of H 2 0 
throughout an asteroidal interior and the eventual loss of volatiles to space. In this 
appendix, a simple dimensional analysis of the degassing from the asteroidal 
surface is given, and a relation is derived for the depth to which free H 2 0 should be 
largely depleted as a function of orbital semimajor axis and regolith permeability. 

Consider the one-dimensional diffusion of a constant-viscosity gas at low 
temperature in a medium of fixed permeability. Gas transport must satisfy the 
equation of state, the continuity equation, and Darcy's law: 


PV = nKT 


d(pu) _ dp 
dz 0t 


KdP 
n dz 


(6.B1) 


where all symbols follow the definitions given in Appendix 6. A. Equations (6.B1) 
may be combined into the single equation 


9P = K 
9t tj 



(6.B2) 


Dimensional analysis of (6.B2) gives, to order of magnitude, 


2 

AP_ K AP 
At Tj Az 


(6.B3) 
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from which follows the time scale t 0 for degassing 


l 0' 


_ T| Az 
KAP 


(6.B4) 


Equation (6.B4) applies to a region without gas sources. However, if the interior 
is in a state of vapor-condensate equilibrium, a special, self-limiting source exists: 
the vapor pressure P vap is maintained as long as condensate is available. Under 
equilibrium conditions, the amount of vapor, expressed in moles, is just n vap = 
v v P vap V//7T, where v v is the volume fraction of voids (~ 1 - v r -Vj). Since asteroidal 
interiors are nearly isothermal over most of solar system history, the vapor pressure 
may be evaluated at the surface temperature, and the pressure drop from the surface 
through the degassed zone to the region of vapor-condensate equilibrium is AP = 
P vap . From (6.B4), then, an amount of H 2 0 ~ n vap is lost on a time scale t 0 . However, 
the total amount of H 2 0 present is n tot = Vj pj V / M, where M is the molar mass. 
Therefore the number of characteristic degassing time scales that can be endured 
before ice is exhausted is 

N = n tot / n vap = Vj pj R T / v v M Pyap (6.B5) 


Dropping the comparable quantities Vj and v v , the total degassing time scale, t = 
N t 0 , is 


t 


p i T| R T Az 


KMP 


vap 


(6.B6) 
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which may be rearranged to give the depth scale Az of degassing over time scale t: 


Az = P (T) 
vap v ' 



(6.B7) 


Substituting (6.A10) for the vapor pressure and (6.A5) for the surface temperature, Az 
may be found as a function of permeability and orbital semimajor axis. Figure 6.B1 
compares the approximate result (6.B7) with a numerical calculation performed 
using the analysis given in Appendix 6.A, with K = 10' 7 m 2 > rj = 1 0 -5 Pa-s, and at a 
distance of 1 A.U. Reasonable agreement is found for the depth of the devolatized 
zone as a function of time. In contrast, the time scales for similar degassing depths 
given by (6.B4) or by the rate of free sublimation [e.g., Zemansky and Dittman, 1981] 
are only a month or less. 
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Table 6.A1. Assumed abundances of heat-producing elements 


Nuclide 

Initial Abundance a 
(mass fraction) 

Decay Constant 

(r‘) 

Initial Heat Production 
(W/kg) 

26A1 

(^Al/^Al) specified 

9.63xl0- 7 

4.21xl0- 3 (26A1/27A1) b 

238u 

2.64xl0* 8 

1.55xl0- 10 

2.48xl0' 12 

235u 

8.60x10-9 

7.91x10-9 

4.98xl0- 12 

232-rh 

5.21x10-8 

4.99xl0* n 

1.36xl0- 12 

40 K 

6.60x10-7 

SASxIO-w 

19.9xl0 12 


a Observed C2 abundances [Mason, 1971] extrapolated to t = 4.60 b.y. before present, 
b Using observed 27 A1 abundance of 1.17 wt% in CM chondrites [Dodd, 1981]. 
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FIGURE CAPTIONS 

Figure 6.1 . Schematic illustration of the alternative models tested in this paper for 
aqueous alteration of carbonaceous chondrites, (a) Alteration occurs 
uniformly throughout the parent body interior. Hydrothermal circu- 
lation may develop at sufficiently large thermal gradients, (b) 
Alteration occurs within a regolith. Under internal heating, water may 
be supplied to the regolith by in situ melting and hydrothermal 
circulation, by venting through fractures that accompany tensile failure 
of the parent body interior at high fluid or gas pressures, or by simple 
vapor diffusion through pre-existing pores and cracks. Impacts can 
also contribute to melting of ice within the regolith. 

Figure 6.2. Profiles of temperature versus radius within a parent body 100 km in 
diameter with 30% ice by volume and an initial 26 AI/ 27 AI ratio of 

l. 6x1 O' 6 . The deep interior is brought just past H 2 0 melting in 10 

m. y. This scenario gives the minimum heating required to produce 
interior aqueous alteration in at the specified size and composition. 

Figure 6.3. Thermal evolution of an object similar to that in Figure 6.2, but with 
initial 26 AI/ 27 ’AI = 2.0x1 O' 6 . The average peak temperature of all rocks 
exposed to liquid water is 25°C, matching the peak temperature of 
CM chondrites inferred from oxygen isotope studies [Clayton and 
Mayeda, 1984]. However, the altered rock fraction is not 
representative of the interior. 
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Figure 6.4. Thermal evolution of an object similar to that in Figure 6.2, but 
matching the Cl peak temperature of 150°C and illustrating potential 
effects of hydrothermal circulation, (a) Initial 26 AI/ 27 AI = 4.3x1 O' 6 ; 
convection suppressed. Altered fraction is representative of the 
interior, (b) Initial 26 AI/ 27 AI = 6.0x1 O' 6 , with vigorous convection 
allowed so that excess heat is removed. The approximation for 
convective heat transfer yields an overestimate of the Nusselt 
number, so the indicated heat production is an upper limit at this size 
and composition. 

Figure 6.5. Summary of allowed initial 26 AI/2 7 AI ratio as a function of parent body 
size, ice:rock composition, and carbonaceous chondrite group 
produced. In each panel, the lower curve gives the minimum 
26 AI/ 27 AI required just to melt interior ice. The next two curves give 
initial 26 AI/ 27 AI required to match peak temperatures of CM and Cl 
chondrites, respectively. The uppermost (dashed) curve gives an 
upper bound to 26 AI abundance in Cl objects under efficient 
hydrothermal circulation; gravitational instability of water does not 
occur for CM models. Long-lived radionuclides alone are sufficient to 
heat D = 1000 object to Cl temperatures. Open square gives maxi- 
mum initial 26 AI/ 27 AI for anhydrous CO chondrites. Comparable 
26 AI/ 27 AI ratios inferred for Cl and CO chondrite parent bodies and 
also for ordinary chondrite parent objects [Miyamoto et at., 1981; 
Hutcheon et ai, 1988] suggest that water may have played an 
important role in their different thermal histories. 
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Figure 6.6. Evolution of a CM object (D = 100 km, ice fraction = 0.2, initial 
26 AI/ 27 AI = 1 .2x1 0* 6 ) explicitly including thermal and compositional 
changes associated with aqueous alteration, (a) Temperature 
profiles. Following complete absorption of the latent heat of fusion of 
H 2 0, temperatures increase sharply due to release of chemical heat 
of reaction. Oxygen isotopic constraints apply only during aqueous 
alteration, so peak temperature is limited instead by inferences from 
carbonaceous chondrite textures, (b) Compositional profile following 
alteration, showing volume fractions of H 2 0 (solid line), anhydrous 
minerals (long dash), and hydrous minerals (short dash). Hydrous: 
anhydrous mineral ratio matches CM proportions, but altered zone is 
confined to the deep interior in this example. 

Figure 6.7. Evolution of an object similar to that in Figure 6.6, but with initial 
26AI/ 27 AI = 5x1 0* 6 . (a) Temperature profiles, showing rapid rise 
through hydration at 273 K to dehydration temperature (623K) and 
maintenance of near-constant value during absorption of the latent 
heat of reaction. High pore pressures due to release of superheated 
water lead to failure of the surrounding rock and venting of H 2 0. (b) 
Final compositional profile. Venting of H 2 0 following dehydration of 
deep interior leaves only anhydrous minerals there. Altered zone is 
restricted to 1-1 1 km depth. 

Figure 6.8. Object with the same size and initial 26 AI/ 27 AI ratio as Figure 6.7, but 
initial H 2 0 is doubled (volume fraction = 0.4). Such an H 2 0 fraction is 
needed to achieve the nearly complete alteration of anhydrous 
minerals observed in Cl chondrites. Heating is insufficient to drive 
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dehydration reaction, so dehydration temperature can act as a barrier 
to further temperature excursions. The final state (not shown) consists 
nearly uniformly of altered minerals, although in reality partial 
dehydration textures would be evident. 

Figure 6.9. Aqueous alteration within the regolith of a large (D=300 km) parent 
body with initial 26 AI /2 7 AI * 5x1 0* 6 , but the pore volume of 25% is 
taken to be undersaturated with only 10% H 2 0. Hydrothermal 
circulation and chemical reactions are not included, (a) 
Temperature profiles. If only regolith is sampled, the high internal 
temperatures are irrelevant, (b) Final H 2 0 distribution. Large internal 
gas pressures cause failure of the surrounding rock, leading to 
venting of 70% of the object's H 2 0. If part of this vented H 2 0 
condenses in the regolith, a favorable environment for aqueous 
alteration would be provided. 

Figure 6.10. Regolith-alteration calculation following that in Figure 6.9, but 
allowing hydrothermal circulation. Vigorous convection holds down 
temperatures in the interior so that only minor venting of H 2 0 occurs. 
Melting also extends to within 1 km of the surface, closer than in 
Figure 6.9. Regolithic water is continuously recharged by hydro- 
thermal circulation. 

Figure 6.11. Regolith-alteration model following that in Figure 6.9, but including 
H 2 0 vapor diffusion, (a) Temperature profiles, (b) H 2 0 distribution. 
Hydrothermal convection must be suppressed or else water 
distribution is homogenized. Upper 10 km has steady supply of vapor 
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over time scales of ~ 1 m.y., providing an additional mechanism for 
introduction of water into the regolith. 

Figure 6.12. Depth of H 2 0 loss over the age of the solar system as a function of 
reference permeability, for different values of solar distance. Effective 
permeability may be enhanced by low-temperature molecular flow 
and depends on pore size, (a) Pore diameter 1 mm, appropriate to 
the sizes of chondrules and clasts in carbonaceous chondrites, (b) 
Pore diameter 1 pm, appropriate to the size of carbonaceous 
chondrite matrix particles. Easier transition to molecular flow occurs 
for smaller pores, allowing greater H 2 0 loss. Degassing at asteroidal 
distances is nevertheless limited to depths of tens of meters. 

Figure 6.B1 . Example of one-dimensional vapor diffusion at 1 A.U. solar distance, 
for a regolith with 10% H 2 0 volume, 25% porosity, and permeability 
10*7 m 2 . Numerical calculation (solid) compares well with simple 
approximation (dash) from dimensional analysis for the depth of 
nearly complete H 2 0 loss. 
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CHAPTER 7: 

CONCLUDING REMARKS 
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In this thesis, spacecraft measurements and laboratory analyses have been 
combined to infer features of the structure and evolution of the lithosphere of Venus 
and the parent bodies of chondritic meteorites. In Part I, limits to both the present 
thickness and the accumulation rate of crust on Venus have been obtained from 
the morphology and distribution of impact craters. Radar altimetry and imaging 
have been used to test hypotheses for crustal divergence. In Part II, models for the 
early geologic histories of both ordinary and carbonaceous chondrite parent 
bodies have been developed in accordance with constraints imposed by chemical, 
isotopic, and textural data. 

The principal conclusions of Part I are 

(1) The dominant mode of lithospheric heat loss on Venus is conduction. 
Volcanism delivers a negligible fraction of the planetary heat loss. Litho- 
spheric recycling cannot be ruled out but can account for no more than 
about one quarter of the planetary heat budget. 

(2) The mean crustal thickness throughout the lowlands and rolling plains is no 
more than 10-20 km. This implies either a low rate of crustal production or 
some form of crustal recycling. 

(3) The question of horizontal crustal divergence in Aphrodite Terra is 
unresolved. There is no quantitative evidence to support suggestions that 
individual topographic elements hundreds of kilometers in size are 
bilaterally symmetric on a regional scale or that the topography is consistent 
with a diverging and cooling thermal boundary layer. These characteristics 
cannot be completely excluded, but neither can they be invoked to support 
a rigid-plate divergence model. Detailed mapping by the Magellan 
spacecraft will resolve the geologic features needed to discriminate 
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between the vertical-tectonics, plate-divergence, and distributed- 
deformation models. 

From these considerations, we may pose two alternative scenarios for the 
global tectonics of Venus: 

(a) Hot, buoyant lithosphere is unable to subduct and is therefore unable to 
drive significant horizontal motion and lithospheric divergence. Crustal 
production, associated with pressure-release melting at sites of lithospheric 
divergence, is inhibited. The crust has accumulated slowly over geologic 
time. 

(b) Crust and lithosphere are recycled, but at a lower rate than on Earth. The 
specific recycling mechanism is presently unknown. 

Is Venus an active planet, truly Earth's twin? Venus must be tectonically active 
because the correlation of long-wavelength topography and gravity points to a 
strong interaction of mantle convection with the surface, and because, in the 
absence of dynamic support, long wavelength topographic features would be 
removed by viscous relaxation within a few hundred million years. The extent of 
igneous activity depends on the mechanism and rate of crustal recycling under 
scenario (a) or (b). 

The principal conclusions of Part II are 

(1) Collisional fragmentation and gravitational reassembly of meteorite parent 
bodies during metamorphism can account for uncorrelated peak temper- 
atures and cooling rates among ordinary chondrites. This model may be 
preferred over the alternative metamorphosed-planetesimal model because 
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there are fewer thermal and collisional restrictions. Moreover, fragmen- 
tation and reassembly provides a single mechanism to explain both 
narrowly-distributed and widely-scattered cooling rates among breccia 
clasts, by collisions during or after metamorphism, respectively. 

(2) Incorporation of ice into carbonaceous chondrite parent bodies can strongly 
influence the thermal evolution of these objects. If aqueous alteration 
occurred throughout carbonaceous chondrite parent body interiors, water 
may have acted as a thermal buffer, allowing aqueous alteration to replace 
the thermal recrystallization observed for ordinary chondrites. If aqueous 
alteration occurred within a surficial regolith, then water must be supplied 
there by direct melting of local ice (which may then undergo hydrothermal 
circulation), by venting of liquid or vapor along fractures caused by failure 
under high pore pressures, or by vapor diffusion through existing pores and 
cracks. Retention of H 2 0 on parent bodies is not limited by sublimation or 
by collisional comminution, but by shock vaporization. Large C-type 
asteroids, if representative of carbonaceous chondrite parent bodies, might 
have significant quantities of ice in their interiors. 

These studies have addressed aspects of the different geological evolutions of 
the largest and smallest terrestrial objects. Study of the full range of geological 
diversity in the solar system offers the best promise of understanding the origin and 
evolution of the Earth and the future potential for human exploration and 
development. 



Summer grass: 

of stalwart warriors' splendid dreams 
the aftermath. 


Basho, on his visit to an old battlefield. 



